hanalyze-cli (empty) → 0.2.0.1
raw patch · 4 files changed
+4159/−0 lines, 4 filesdep +basedep +containersdep +dataframe-core
Dependencies added: base, containers, dataframe-core, dataframe-csv, dataframe-operations, filepath, hanalyze, hmatrix, hvega, mwc-random, text, time, vector
Files
- README.ja.md +111/−0
- README.md +112/−0
- app/Main.hs +3883/−0
- hanalyze-cli.cabal +53/−0
+ README.ja.md view
@@ -0,0 +1,111 @@+# hanalyze-cli++[`hanalyze`](../README.ja.md) のコマンドラインフロントエンド。+実行ファイル **`hanalyze`** 1 本だけを提供する package で、 library は持たない。++Phase 106 で umbrella package から分離した。 目的は 2 つ:++- library の作業をするたびに CLI の compile + link が走らないようにする+- **umbrella の公開 API だけに依存する**構成にして、 internal module へ手を+ 伸ばしていないことを構造で保証する (`build-depends` は `hanalyze` の+ 1 本のみ)++```bash+cabal build hanalyze # 既定の cabal.project に含まれる+cabal run hanalyze -- --help+```++## サブコマンド (全 15 個)++### 回帰++| コマンド | 内容 |+|---|---|+| `regress` | 古典 / ベイズ回帰 (LM / GLM / GLMM / GP / HBM)。 **既定のコマンド**で、 サブコマンド名を省くとこれになる |+| `ridge` | 罰則付き回帰 (Ridge / Lasso / Elastic Net) |+| `kernel` | カーネル回帰 / RFF 近似 |+| `spline` | B-spline / 自然三次スプライン回帰 |+| `quantile` | 分位点回帰 (τ 分位・MM-IRLS) |+| `gam` | 一般化加法モデル (加法 B-spline + Ridge) |+| `rf` | ランダムフォレスト回帰 (CART + bagging + 特徴量サブセット) |+| `multireg` | 多出力回帰 (wide CSV・linear / kernel-rbf) |++### データの下見・可視化++| コマンド | 内容 |+|---|---|+| `info` | 列ごとの型と基本統計量を表示 |+| `hist` | ヒストグラム (理論分布の重ね描きも可) |++### 実験計画法++| コマンド | 内容 |+|---|---|+| `doe` | 直交表 (L_n) の生成 |+| `taguchi` | 田口メソッド (SN 比 + 要因効果 + 内側 / 外側配置) |++### データ整形++| コマンド | 内容 |+|---|---|+| `clean` | 列ごとのクリーニング規則を適用 (`StripUnits` / `ParseCurrency` / `ParseDecimalEU` …) |+| `melt` | wide → long 変換 |+| `regrid` | 歯抜けの long 形式データ `[id, z, y]` を共通 grid へ揃える |++> ⚠ `clean` / `melt` / `regrid` は実装済みだが `hanalyze --help` の一覧には+> 載っていない。 各コマンドを引数なしで実行すれば usage が出る+> (例: `hanalyze melt`)。++## 使い方++```bash+# 列ごとの型・基本統計を見る+hanalyze info data.csv++# 単回帰 (サブコマンド省略 = regress)+hanalyze data.csv x y++# 多項式の次数を指定して 90% 信頼区間つき+hanalyze data.tsv "x1 x2" y LM --degree -1 2 -2 3 --ci 0.90++# ポアソン回帰 (log link)+hanalyze data.csv x y GLM -d poisson -l log++# 変量効果つき (LM + --group → LME、 GLM + --group → GLMM)+hanalyze data.csv x y LM --group school++# HTML レポートを書き出す+hanalyze data.csv x y --report report.html --waic+```++`regress` の主なオプション:++| オプション | 内容 |+|---|---|+| `-d, --dist DIST` | 分布 `gaussian` / `binomial` / `poisson` (既定 `gaussian`) |+| `-l, --link LINK` | リンク関数 `identity` / `log` / `logit` / `sqrt` (既定 = 正準リンク) |+| `--degree SPEC` | 多項式次数の指定 (既定 `1`) |+| `--ci [LEVEL]` | 信頼区間 (既定 `0.95`) |+| `--pi [LEVEL]` | 予測区間 (Gaussian のみ・既定 `0.95`) |+| `--group COL` | グループ列 → LME / GLMM |+| `--format FORMAT` | 出力形式 `html` / `png` / `svg` (既定 `html`) |+| `--report [FILE]` | HTML 解析レポートを生成 (既定 `report.html`) |+| `--waic` | WAIC と LOO-CV をレポートに載せる (`--report` が前提) |++各サブコマンド固有のオプションは、 引数なしで実行すると表示される+(例: `hanalyze ridge`)。++## 関連 docs++CLI 専用のページはまだ無く、 各機能の doc の中で `hanalyze <sub>` の実行例が+示されている。++- 回帰の入口: [docs/regression/01-lm.ja.md](../docs/regression/01-lm.ja.md)+- `clean` の規則一覧: [docs/io/01-dirty-data.ja.md](../docs/io/01-dirty-data.ja.md)+- `melt` / `regrid`: [io/02-reshape.ja.md](../docs/io/02-reshape.ja.md) /+ [io/03-regrid.ja.md](../docs/io/03-regrid.ja.md)+- `doe` / `taguchi`: [doe/01-doe.ja.md](../docs/doe/01-doe.ja.md) /+ [doe/02-orthogonal-taguchi.ja.md](../docs/doe/02-orthogonal-taguchi.ja.md)+- `--report` が出す HTML: [visualization/02-report-builder.ja.md](../docs/visualization/02-report-builder.ja.md)++← [repository README](../README.ja.md)
+ README.md view
@@ -0,0 +1,112 @@+# hanalyze-cli++The command-line front end of [`hanalyze`](../README.md). This package+ships a single executable, **`hanalyze`**, and no library.++It was split out of the umbrella package in Phase 106 for two reasons:++- so that working on the library does not trigger a CLI compile + link every+ time, and+- so that the CLI **depends only on the umbrella's public API** — the fact+ that it never reaches into internal modules is guaranteed structurally+ (`build-depends` names only `hanalyze`).++```bash+cabal build hanalyze # included in the default cabal.project+cabal run hanalyze -- --help+```++## Subcommands (15 in total)++### Regression++| Command | What it does |+|---|---|+| `regress` | Classical / Bayesian regression (LM / GLM / GLMM / GP / HBM). **The default** — omitting the subcommand name lands here |+| `ridge` | Regularized regression (Ridge / Lasso / Elastic Net) |+| `kernel` | Kernel regression / RFF approximation |+| `spline` | B-spline / natural cubic regression |+| `quantile` | Quantile regression (τ-quantile, MM-IRLS) |+| `gam` | Generalized additive model (additive B-splines + Ridge) |+| `rf` | Random forest regression (CART + bagging + feature subsets) |+| `multireg` | Multi-output regression (wide CSV; linear / kernel-rbf) |++### Inspecting and plotting data++| Command | What it does |+|---|---|+| `info` | Print per-column type and basic statistics |+| `hist` | Histogram, optionally with a theoretical density overlaid |++### Design of experiments++| Command | What it does |+|---|---|+| `doe` | Generate orthogonal arrays (L_n) |+| `taguchi` | Taguchi method (S/N ratio + factor effects + inner/outer arrays) |++### Reshaping data++| Command | What it does |+|---|---|+| `clean` | Apply per-column cleaning rules (`StripUnits`, `ParseCurrency`, `ParseDecimalEU`, …) |+| `melt` | Wide → long reshaping |+| `regrid` | Align sparse long-form data `[id, z, y]` onto a common grid |++> ⚠ `clean`, `melt` and `regrid` are implemented but missing from the+> `hanalyze --help` listing. Run any of them with no arguments to see its+> usage (e.g. `hanalyze melt`).++## Usage++```bash+# inspect column types and basic statistics+hanalyze info data.csv++# simple regression (subcommand omitted = regress)+hanalyze data.csv x y++# polynomial degrees with a 90% confidence interval+hanalyze data.tsv "x1 x2" y LM --degree -1 2 -2 3 --ci 0.90++# Poisson regression with a log link+hanalyze data.csv x y GLM -d poisson -l log++# random effects (LM + --group → LME, GLM + --group → GLMM)+hanalyze data.csv x y LM --group school++# write an HTML report+hanalyze data.csv x y --report report.html --waic+```++Main options of `regress`:++| Option | Meaning |+|---|---|+| `-d, --dist DIST` | Distribution: `gaussian` / `binomial` / `poisson` (default `gaussian`) |+| `-l, --link LINK` | Link: `identity` / `log` / `logit` / `sqrt` (default: canonical) |+| `--degree SPEC` | Polynomial degree specification (default `1`) |+| `--ci [LEVEL]` | Confidence interval (default `0.95`) |+| `--pi [LEVEL]` | Prediction interval (Gaussian only; default `0.95`) |+| `--group COL` | Grouping column → LME / GLMM |+| `--format FORMAT` | Output format: `html` / `png` / `svg` (default `html`) |+| `--report [FILE]` | Generate an HTML analysis report (default `report.html`) |+| `--waic` | Add WAIC and LOO-CV to the report (requires `--report`) |++Options specific to the other subcommands are printed when you run them with+no arguments (e.g. `hanalyze ridge`).++## Related docs++There is no CLI-specific page yet; `hanalyze <sub>` examples appear inside the+per-feature docs.++- Regression overview: [docs/regression/01-lm.md](../docs/regression/01-lm.md)+- `clean` rules: [docs/io/01-dirty-data.md](../docs/io/01-dirty-data.md)+- `melt` / `regrid`: [io/02-reshape.md](../docs/io/02-reshape.md) /+ [io/03-regrid.md](../docs/io/03-regrid.md)+- `doe` / `taguchi`: [doe/01-doe.md](../docs/doe/01-doe.md) /+ [doe/02-orthogonal-taguchi.md](../docs/doe/02-orthogonal-taguchi.md)+- The HTML that `--report` produces: [visualization/02-report-builder.md](../docs/visualization/02-report-builder.md)++← [repository README](../README.md)
+ app/Main.hs view
@@ -0,0 +1,3883 @@+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE RankNTypes #-}+module Main where++import Hanalyze.DataIO.CSV (loadAutoSafeWith, LoadOpts (..), defaultLoadOpts)+import qualified Hanalyze.DataIO.Log as Log+import qualified Hanalyze.DataIO.Clean as Clean+import qualified Hanalyze.Stat.Standardize as Std+import qualified Hanalyze.Stat.NumberFormat as NF+import Data.Time.Clock (getCurrentTime, diffUTCTime, UTCTime)+import qualified Hanalyze.DataIO.Preprocess as Pp+import qualified Hanalyze.Stat.Interpolate as Interp+import qualified Hanalyze.Stat.AdaptiveGrid as AG+import Text.Read (readMaybe)+import qualified DataFrame.Internal.DataFrame as DX+import qualified DataFrame.Operations.Core as DX+import qualified DataFrame.IO.CSV as DX+import qualified DataFrame.Internal.Column as DXC+import qualified DataFrame.Internal.DataFrame as DXD+import Hanalyze.DataIO.Convert (getDoubleVec, getTextVec, getMaybeTextVec)+import Hanalyze.Model.Core (Band (..), FitResult, rSquared1, coeffList, fittedList, residualsV)+import qualified Hanalyze.Model.Core as Core+import Hanalyze.Model.GLM (Family (..), parseFamily, LinkFn (..), parseLink, canonicalLink,+ fitGLMWithSmooth, fitGLMFull)+import Hanalyze.Model.GLMM (GLMMResult (..), fitLMEDataFrame, fitGLMMDataFrame)+import Hanalyze.Model.LM (SmoothFit (..), multiPolyDesignMatrix)+import Hanalyze.Stat.Distribution (Distribution, parseDistribution)+import Hanalyze.Viz.Core (defaultConfig, openInBrowser, OutputFormat (..), parseFormat)+import Hanalyze.Viz.Scatter (scatterWithSmoothFile, scatterMultiYFile, scatterPlotFile,+ scatterWithGroupsFile, predictedVsActualFile,+ predictedVsActual, scatterWithGroups)+import Hanalyze.Viz.Histogram (histogramPlotFile, histogramWithDensityFile)+import Hanalyze.Viz.AnalysisReport (AnalysisReportConfig (..), ModelFit (..), NamedPlot (..),+ SmoothData (..), GPKernelFit (..), GPFitSummary (..), FitSummary (..),+ GLMMSummary (..), HBMRegSummary (..),+ mkFitSummary, mkGLMMSummary,+ writeAnalysisReport, writeAnalysisReportPlots)+import qualified Hanalyze.Design.Orthogonal as OA+import qualified Hanalyze.Design.Taguchi as TG+import qualified Hanalyze.Viz.Taguchi as VTG+import qualified Hanalyze.Viz.ReportBuilder as RB+import qualified Hanalyze.Viz.ReportInstances as RI+import qualified Hanalyze.Viz.ModelGraph+import qualified Graphics.Vega.VegaLite as VL+import Graphics.Vega.VegaLite (VegaLite, VLProperty, VLSpec)+import qualified Hanalyze.Model.KernelRegression as Kern+import qualified Hanalyze.Model.MultiLM as MLM+import qualified Hanalyze.Model.Regularized as Reg+import qualified Hanalyze.Model.GAM as GAM+import qualified Hanalyze.Model.Quantile as QR+import qualified Hanalyze.Model.RandomForest as RF+import qualified Hanalyze.Model.RFF as RFF+import qualified Hanalyze.Model.Spline as Spl+import Hanalyze.Model.LM (SmoothFit (..))+import qualified Hanalyze.Model.HBM as HBMod+import qualified Hanalyze.MCMC.NUTS as HBMnuts+import qualified Hanalyze.MCMC.Core as MCMCcore+import qualified Data.Map.Strict as Map+import Hanalyze.Viz.MCMC (mcmcDiagnostics, autocorrPlot)+import Hanalyze.Viz.Core (PlotConfig (..))+import Hanalyze.Model.GP (Kernel (..), GPModel (..), GPParams, GPPredData,+ GPResult, gpMean,+ initParamsFromData, optimizeGP, fitGP, logMarginalLikelihood,+ gpPredData)++import Hanalyze.Stat.ModelSelect (lmPosteriorLogLiks, glmPosteriorLogLiks,+ lmePosteriorLogLiks, waic, loo,+ WAICResult (..), LOOResult (..))++import Control.Monad (when)+import Data.Char (isDigit)+import Data.List (intercalate, sort)+import qualified Data.Set as Set+import System.FilePath (dropExtension)+import qualified Data.Text as T+import qualified Data.Text.IO as TIO+import qualified Data.Vector as V+import qualified Numeric.LinearAlgebra as LA+import System.Environment (getArgs)+import System.IO (hPutStrLn, stderr)+import System.Random.MWC (createSystemRandom)+import Text.Printf (printf)++-- ---------------------------------------------------------------------------+-- CLI types+-- ---------------------------------------------------------------------------++data ModelType = LM | GLM | NoReg | GP | HBM deriving (Show, Eq)++data DegreeSpec+ = AllDegree Int+ | PerDegree [(Int, Int)]+ deriving (Show)++data Config = Config+ { cfgFile :: FilePath+ , cfgXCols :: [T.Text]+ , cfgYCols :: [T.Text] -- one or more y columns+ , cfgModel :: ModelType+ , cfgDist :: Family+ , cfgLink :: LinkFn+ , cfgDegree :: DegreeSpec+ , cfgBand :: Band+ , cfgFormat :: OutputFormat+ , cfgGroup :: Maybe T.Text -- grouping column → LME / GLMM+ , cfgHistMode :: Bool -- --hist: draw histogram of x column+ , cfgFitDist :: Maybe Distribution -- --fit DIST PARAMS+ , cfgReport :: Maybe FilePath -- --report [FILE]: generate HTML report+ , cfgWAIC :: Bool -- --waic: compute WAIC/LOO-CV+ , cfgLoadOpts :: LoadOpts -- --no-header / --skip / --comment / --strict+ } deriving (Show)++-- ---------------------------------------------------------------------------+-- Argument parsing+-- ---------------------------------------------------------------------------++usageMsg :: String+usageMsg = unlines+ [ "Usage: hanalyze <file> <xcols> <ycols> [LM|GLM|NoReg|GP|HBM] [options]"+ , ""+ , " <file> CSV/TSV/SSV file (auto-detected from extension)"+ , " <xcols> x column name(s); quote multiple: \"x1 x2\""+ , " <ycols> y column name(s); quote multiple: \"y1 y2\" (multi-y → scatter only)"+ , " LM|GLM|NoReg|GP|HBM model type (default: LM)"+ , " GP: Gaussian Process regression (single x/y only); compares RBF, Matérn5/2, Periodic"+ , " HBM: Bayesian linear regression via NUTS (single x/y only); --report で AnalysisReport 生成"+ , ""+ , "Options:"+ , " -d, --dist DIST distribution: gaussian|binomial|poisson (default: gaussian)"+ , " -l, --link LINK link function: identity|log|logit|sqrt (default: canonical)"+ , " --degree SPEC degree specification (default: 1)"+ , " --ci [LEVEL] show confidence interval (default level: 0.95)"+ , " --pi [LEVEL] show prediction interval (Gaussian only; default level: 0.95)"+ , " --format FORMAT output format: html|png|svg (default: html)"+ , " --group COL grouping column → LM+group: LME, GLM+group: GLMM"+ , " --report [FILE] generate HTML analysis report (default: report.html)"+ , " --format png|svg と組み合わせるとプロット部分を画像にも出力"+ , " --waic compute WAIC and LOO-CV and show in report (requires --report)"+ , ""+ , "Degree specification:"+ , " N all columns get degree N"+ , " -i1 N1 [-i2 N2…] column at 1-based position i1 gets degree N1; others: 1"+ , ""+ , "Examples:"+ , " hanalyze data.csv x y"+ , " hanalyze data.tsv \"x1 x2\" y LM --degree -1 2 -2 3 --ci 0.90"+ , " hanalyze data.csv x y GLM -d poisson -l log"+ , " hanalyze data.csv x y LM --group school"+ , " hanalyze data.csv x y GLM -d binomial -l logit --group hospital"+ , " hanalyze data.csv x \"y1 y2\" NoReg"+ ]++parseArgs :: [String] -> Either String Config+parseArgs args0 =+ let (lopts, args) = parseLoadOpts args0+ in case args of+ (file : xColsStr : yColsStr : rest) -> do+ let xCols = map T.pack (words xColsStr)+ yCols = map T.pack (words yColsStr)+ if null xCols+ then Left "Error: xcols must not be empty"+ else if null yCols+ then Left "Error: ycols must not be empty"+ else do+ (model, rest1) <- parseModelType rest+ (mDist, mLink, degSpec, band, fmt, mGrp, hist, mFit, mRpt, waicF, rest2) <- parseOptions rest1+ if not (null rest2)+ then Left ("Unexpected argument(s): " ++ unwords rest2)+ else do+ let dist = maybe Gaussian id mDist+ lnk = maybe (canonicalLink dist) id mLink+ Right Config+ { cfgFile = file+ , cfgXCols = xCols+ , cfgYCols = yCols+ , cfgModel = model+ , cfgDist = dist+ , cfgLink = lnk+ , cfgDegree = degSpec+ , cfgBand = band+ , cfgFormat = fmt+ , cfgGroup = mGrp+ , cfgHistMode = hist+ , cfgFitDist = mFit+ , cfgReport = mRpt+ , cfgWAIC = waicF+ , cfgLoadOpts = lopts+ }+ _ -> Left usageMsg++parseModelType :: [String] -> Either String (ModelType, [String])+parseModelType ("LM" : rest) = Right (LM, rest)+parseModelType ("GLM" : rest) = Right (GLM, rest)+parseModelType ("NoReg" : rest) = Right (NoReg, rest)+parseModelType ("GP" : rest) = Right (GP, rest)+parseModelType ("HBM" : rest) = Right (HBM, rest)+parseModelType rest = Right (LM, rest)++parseOptions :: [String]+ -> Either String (Maybe Family, Maybe LinkFn, DegreeSpec, Band, OutputFormat,+ Maybe T.Text, Bool, Maybe Distribution, Maybe FilePath, Bool, [String])+parseOptions = go Nothing Nothing (AllDegree 1) NoBand HTML Nothing False Nothing Nothing False+ where+ go mDist mLink deg band fmt mGrp hist mFit mRpt waicF [] =+ Right (mDist, mLink, deg, band, fmt, mGrp, hist, mFit, mRpt, waicF, [])++ go mDist mLink deg band fmt mGrp hist mFit mRpt waicF (flag : rest)+ | flag `elem` ["-d", "--dist"] = case rest of+ (v:rest') -> do fam <- parseFamily v+ go (Just fam) mLink deg band fmt mGrp hist mFit mRpt waicF rest'+ [] -> Left "Error: -d/--dist requires an argument"++ | flag `elem` ["-l", "--link"] = case rest of+ (v:rest') -> do lnk <- parseLink v+ go mDist (Just lnk) deg band fmt mGrp hist mFit mRpt waicF rest'+ [] -> Left "Error: -l/--link requires an argument"++ | flag == "--degree" = do+ let (degTokens, remaining) = span isDegreeToken rest+ if null degTokens+ then Left "--degree requires a specification (e.g., 2 or -1 2 -2 3)"+ else do degSpec <- parseDegreeSpec degTokens+ go mDist mLink degSpec band fmt mGrp hist mFit mRpt waicF remaining++ | flag == "--ci" =+ let (level, rest') = consumeLevel 0.95 rest+ in go mDist mLink deg (CI level) fmt mGrp hist mFit mRpt waicF rest'++ | flag == "--pi" =+ let (level, rest') = consumeLevel 0.95 rest+ in go mDist mLink deg (PI level) fmt mGrp hist mFit mRpt waicF rest'++ | flag `elem` ["-f", "--format"] = case rest of+ (v:rest') -> do f <- parseFormat v+ go mDist mLink deg band f mGrp hist mFit mRpt waicF rest'+ [] -> Left "Error: -f/--format requires an argument"++ | flag == "--group" = case rest of+ (v:rest') -> go mDist mLink deg band fmt (Just (T.pack v)) hist mFit mRpt waicF rest'+ [] -> Left "Error: --group requires a column name"++ | flag == "--hist" =+ go mDist mLink deg band fmt mGrp True mFit mRpt waicF rest++ | flag == "--fit" = case rest of+ (name:rest') ->+ let (paramStrs, rest'') = span isNumericToken rest'+ params = map read paramStrs :: [Double]+ in case parseDistribution name params of+ Left err -> Left ("--fit: " ++ err)+ Right d -> go mDist mLink deg band fmt mGrp hist (Just d) mRpt waicF rest''+ [] -> Left "--fit requires a distribution name (e.g. --fit normal 0 1)"++ | flag == "--report" = case rest of+ (v:rest') | not (null v) && head v /= '-' ->+ go mDist mLink deg band fmt mGrp hist mFit (Just v) waicF rest'+ _ -> go mDist mLink deg band fmt mGrp hist mFit (Just "report.html") waicF rest++ | flag == "--waic" =+ go mDist mLink deg band fmt mGrp hist mFit mRpt True rest++ | otherwise = Right (mDist, mLink, deg, band, fmt, mGrp, hist, mFit, mRpt, waicF, flag : rest)++isNumericToken :: String -> Bool+isNumericToken s = case (reads s :: [(Double, String)]) of+ [(_, "")] -> True+ _ -> False++-- Consume an optional level (0 < v < 1) after a band flag.+consumeLevel :: Double -> [String] -> (Double, [String])+consumeLevel _ (t:ts) | isLevelToken t = (read t, ts)+consumeLevel def rest = (def, rest)++isLevelToken :: String -> Bool+isLevelToken s = case (reads s :: [(Double, String)]) of+ [(v, "")] | v > 0, v < 1 -> True+ _ -> False++isDegreeToken :: String -> Bool+isDegreeToken ('-' : ds) = not (null ds) && all isDigit ds+isDegreeToken s = not (null s) && all isDigit s++parseDegreeSpec :: [String] -> Either String DegreeSpec+parseDegreeSpec [n] =+ case (reads n :: [(Int, String)]) of+ [(v, "")] | v >= 0 -> Right (AllDegree v)+ _ -> Left ("Invalid degree: " ++ n)+parseDegreeSpec tokens = fmap PerDegree (parsePairs tokens)+ where+ parsePairs [] = Right []+ parsePairs (pos : deg : rest) =+ case (reads pos :: [(Int,String)], reads deg :: [(Int,String)]) of+ ([(p,"")], [(d,"")]) | p < 0, d >= 0 ->+ fmap ((abs p, d) :) (parsePairs rest)+ _ -> Left ("Invalid degree pair near: " ++ pos ++ " " ++ deg)+ parsePairs [t] = Left ("Odd number of tokens in --degree near: " ++ t)++applyDegreeSpec :: DegreeSpec -> [T.Text] -> [(T.Text, Int)]+applyDegreeSpec (AllDegree d) cols = [(c, d) | c <- cols]+applyDegreeSpec (PerDegree ps) cols =+ [ (c, maybe 1 id (lookup i ps)) | (i, c) <- zip [1..] cols ]++-- | --format PNG/SVG が指定されていれば、AnalysisReport のプロットを+-- 個別画像として書き出す (HTML 本体に加えて補助出力)。+maybeExportReportPlots :: Config -> FilePath -> [NamedPlot] -> IO ()+maybeExportReportPlots cfg htmlPath plots =+ case cfgFormat cfg of+ HTML -> return ()+ fmt -> do+ let prefix = dropExtension htmlPath+ paths <- writeAnalysisReportPlots prefix fmt plots+ mapM_ (\p -> putStrLn $ "Plot image: " ++ p) paths++-- ---------------------------------------------------------------------------+-- CLI report builders (Phase 2: regress --report → ReportBuilder 経路)+-- ---------------------------------------------------------------------------++-- | NamedPlot を ReportSection に変換 (タイトル付き secVega)。+namedPlotsToSecs :: [NamedPlot] -> [RB.ReportSection]+namedPlotsToSecs nps =+ [ RB.secVega title vega | NamedPlot _ title vega <- nps ]++-- | WAIC/LOO 結果 (オプション) を 1 セクションに整形。+waicSection :: Maybe (WAICResult, LOOResult) -> [RB.ReportSection]+waicSection Nothing = []+waicSection (Just (w, l)) =+ [ RB.secKeyValue "モデル選択 (WAIC / LOO-CV)"+ [ ("WAIC", T.pack (printf "%.2f" (waicValue w)))+ , ("LOO", T.pack (printf "%.2f" (looValue l)))+ , ("p_WAIC", T.pack (printf "%.2f" (waicPwaic w)))+ , ("k\x0302 > 0.7", T.pack (show (looKHatBad l) ++ " 件"))+ ]+ ]++-- | 残差の (σ_hat, RMSE, max|r|)。p は推定パラメータ数 (intercept 含む)。+cliResidStats :: [Double] -> Int -> (Double, Double, Double)+cliResidStats resid p =+ let n = length resid+ sumSq = sum [ r * r | r <- resid ]+ sH = sqrt (sumSq / fromIntegral (max 1 (n - p)))+ rmse = sqrt (sumSq / fromIntegral (max 1 n))+ mAbs = maximum (0 : map abs resid)+ in (sH, rmse, mAbs)++-- | LM / GLM 用 CLI レポートセクション群。多項式次数と WAIC/LOO に対応。+cliRegressSections+ :: Config -> DXD.DataFrame -> Family -> LinkFn+ -> [(T.Text, Int)] -> FitResult -> Maybe SmoothFit+ -> Maybe (WAICResult, LOOResult)+ -> [NamedPlot]+ -> [RB.ReportSection]+cliRegressSections cfg df dist lnk colDegs res mSmooth mModelSel pvsaPlots =+ let xCols = cfgXCols cfg+ yCol = case cfgYCols cfg of (y:_) -> y; _ -> "y"+ beta = coeffList res+ coefLbls = map T.pack (multiCoeffLabels colDegs)+ coeffs = zip coefLbls beta+ fitted = fittedList res+ resid = LA.toList (residualsV res)+ p = length beta+ (sigmaH, rmse, maxAbs) = cliResidStats resid p+ r2 = rSquared1 res+ r2Lbl = T.pack (r2Label dist)+ isLM = dist == Gaussian+ isPoly = any (\(_, d) -> d > 1) colDegs+ modelType+ | isLM = if isPoly then "LM (polynomial)" else "LM"+ | otherwise = "GLM(" <> T.pack (show dist) <> ")"++ formulaTex+ | isLM = "$" <> yCol <> "_i = "+ <> T.intercalate " + "+ ("\\beta_0" :+ [ "\\beta_" <> T.pack (show (i :: Int)) <> " " <> trm+ | (i, trm) <- zip [1 ..] (polyTerms colDegs) ])+ <> " + \\varepsilon_i$<br>"+ <> "$\\varepsilon_i \\sim \\text{Normal}(0, \\sigma^2)$"+ | otherwise =+ "$g(\\mu_i) = "+ <> T.intercalate " + "+ ("\\beta_0" :+ [ "\\beta_" <> T.pack (show (i :: Int)) <> " " <> trm+ | (i, trm) <- zip [1 ..] (polyTerms colDegs) ])+ <> "$<br>"+ <> "$" <> yCol <> "_i \\sim \\text{" <> T.pack (show dist) <> "}(\\mu_i)$"++ smoothC = case mSmooth of+ Just sf -> RB.SmoothCurve (sfX sf) (sfFit sf) (sfLower sf) (sfUpper sf)+ Nothing -> RB.SmoothCurve [] [] [] []++ scatterCard = case (xCols, mSmooth) of+ ([xc], Just _) -> case (getDoubleVec xc df, getDoubleVec yCol df) of+ (Just xv, Just yv) ->+ [ RB.secCard "散布図 + 回帰線"+ [ RB.secFitScatter xc yCol (V.toList xv) (V.toList yv)+ (Just smoothC) ] ]+ _ -> []+ _ -> []++ -- 対話的予測: 多項式拡張の場合は係数数と x 列数が合わないので省略。+ interactiveSecs = case (isPoly, traverse (`getDoubleVec` df) xCols, getDoubleVec yCol df) of+ (False, Just xVs, Just yV) | not (null xVs) ->+ let xRows = [ [ xv V.! i | xv <- xVs ]+ | i <- [0 .. V.length yV - 1] ]+ mkSlider xv =+ let lo = V.minimum xv+ hi = V.maximum xv+ ext = (hi - lo) * 0.5+ in (lo - ext, (lo + hi) / 2, hi + ext)+ im = RB.InteractiveModel+ { RB.imXCols = xCols+ , RB.imYCol = yCol+ , RB.imXValues = xRows+ , RB.imYValues = V.toList yV+ , RB.imIntercept = head beta+ , RB.imBetas = drop 1 beta+ , RB.imLink = T.pack (linkLabelLower lnk)+ , RB.imSlider = map mkSlider xVs+ , RB.imCISigma = if isLM then Just sigmaH else Nothing+ }+ in [RB.secInteractiveMulti "対話的予測" im]+ _ -> []++ statRow =+ RB.secStatRow+ [ (r2Lbl, T.pack (printf "%.4f" r2))+ , ("方法", if isLM then "OLS (QR)" else "IRLS")+ , ("σ_hat", T.pack (printf "%.4f" sigmaH))+ , ("RMSE", T.pack (printf "%.4f" rmse))+ , ("最大絶対残差", T.pack (printf "%.4f" maxAbs))+ ]++ resultSec =+ RB.secCollapsible "<span class=\"sec-icon\">📈</span> 回帰結果" True+ ([ statRow+ , RB.secCard "係数"+ [RB.secCoefficients coeffs (Just (r2Lbl, r2))]+ ]+ ++ scatterCard+ ++ [RB.secCard "残差プロット" [RB.secResiduals fitted resid]])++ modelSec+ | isLM = RB.secModelOverview modelType formulaTex Nothing+ | otherwise = RB.secModelOverviewLink modelType formulaTex+ (T.pack (linkLabelLower lnk)) Nothing++ extraPlotSecs = namedPlotsToSecs pvsaPlots++ in [ RB.secDataOverview df xCols yCol+ , modelSec+ , resultSec+ ] ++ interactiveSecs ++ extraPlotSecs ++ waicSection mModelSel++-- | colDegs を polynomial 項の文字列に展開: [(x, 2), (z, 1)] → ["x", "x^2", "z"]+polyTerms :: [(T.Text, Int)] -> [T.Text]+polyTerms = concatMap (\(c, d) ->+ [ if k == 1 then c else c <> "^" <> T.pack (show k) | k <- [1 .. d] ])++-- | リンク関数を JS 側のリンク名に対応させる (identity / log / logit / sqrt)。+linkLabelLower :: LinkFn -> String+linkLabelLower Identity = "identity"+linkLabelLower Log = "log"+linkLabelLower Logit = "logit"+linkLabelLower Sqrt = "sqrt"++-- | GLMM (LME) 用 CLI レポートセクション群。+cliMixedSections+ :: Config -> DXD.DataFrame -> Family -> LinkFn+ -> [(T.Text, Int)] -> T.Text -> GLMMResult -> Maybe (WAICResult, LOOResult)+ -> [NamedPlot]+ -> [RB.ReportSection]+cliMixedSections cfg df dist lnk colDegs grpCol gr mModelSel extraPlots =+ let xCols = cfgXCols cfg+ yCol = case cfgYCols cfg of (y:_) -> y; _ -> "y"+ base = RB.toReport (RB.defaultReportConfig "") df xCols yCol+ (RI.GLMMReport gr dist lnk grpCol)+ colDegInfo =+ [ RB.secKeyValue "Polynomial degrees"+ [ (c, T.pack (show d)) | (c, d) <- colDegs, d > 1 ]+ | any (\(_, d) -> d > 1) colDegs+ ]+ in base ++ colDegInfo ++ namedPlotsToSecs extraPlots ++ waicSection mModelSel++-- | GP 用 CLI レポートセクション群。マルチカーネル比較対応。+-- 呼び出し側で予測グリッド X (`gridX`) を渡す。+cliGPSections+ :: T.Text -> T.Text -> DXD.DataFrame -> [Double] -> [Double]+ -> [Double] -- ^ 予測グリッド X+ -> [GPKernelFit]+ -> [RB.ReportSection]+cliGPSections xCol yCol df xs ys gridX kfits =+ let bestK = case kfits of (k:_) -> Just k; _ -> Nothing+ mainSec = case bestK of+ Just kf ->+ RB.toReport (RB.defaultReportConfig "") df [xCol] yCol+ (RI.GPReport (gkKernel kf) (gkParams kf) (gkResult kf)+ gridX xs ys (gkLML kf))+ Nothing -> []+ cmpRows = [ [ gkLabel kf+ , T.pack (printf "%.2f" (gkLML kf)) ]+ | kf <- kfits ]+ cmpSec = case kfits of+ [] -> []+ [_] -> []+ _ -> [RB.secComparisonTable+ "カーネル比較 (LML 降順)"+ ["カーネル", "log p(y|X,θ)"] cmpRows (Just 0)]+ in mainSec ++ cmpSec++-- | HBM 用 CLI レポートセクション群。+cliHBMSections+ :: T.Text -> T.Text -> DXD.DataFrame -> [Double] -> [Double]+ -> MCMCcore.Chain -> Maybe T.Text -> Maybe (WAICResult, LOOResult)+ -> [NamedPlot]+ -> [RB.ReportSection]+cliHBMSections xCol yCol df xs ys chain mGraph mModelSel extraPlots =+ let rep = RI.HBMLinearReport+ { RI.hbmrChain = chain+ , RI.hbmrXs = xs+ , RI.hbmrYs = ys+ , RI.hbmrAlphaName = "alpha"+ , RI.hbmrBetaName = "beta"+ , RI.hbmrSigmaName = "sigma"+ , RI.hbmrGraph = mGraph+ }+ base = RB.toReport (RB.defaultReportConfig "") df [xCol] yCol rep+ _ = xCol -- silence warning+ _ = yCol+ in base ++ namedPlotsToSecs extraPlots ++ waicSection mModelSel++-- ---------------------------------------------------------------------------+-- Main+-- ---------------------------------------------------------------------------++-- ---------------------------------------------------------------------------+-- Subcommand dispatcher (Phase C: hybrid CLI)+--+-- Top-level usage:+-- hanalyze <subcommand> [args...]+-- hanalyze <file> <xcols> <ycols> [LM|GLM|...] [opts] (legacy = regress)+--+-- Implemented: regress, info, hist, help+-- Stubs: ridge, kernel, spline, doe, taguchi+-- ---------------------------------------------------------------------------++helpMsg :: String+helpMsg = unlines+ [ "hanalyze \x2014 general-purpose statistical analysis & visualization toolkit"+ , ""+ , "Usage: hanalyze <subcommand> [args...]"+ , " hanalyze <file> <xcols> <ycols> [LM|GLM|NoReg|GP|HBM] [opts] (legacy = regress)"+ , ""+ , "Subcommands:"+ , " regress Classical/Bayesian regression (LM/GLM/GLMM/GP/HBM) [implemented]"+ , " info Print per-column type and basic statistics [implemented]"+ , " hist Plot a histogram (optionally with theoretical density) [implemented]"+ , " ridge Regularized regression (Ridge/Lasso/Elastic Net) [implemented]"+ , " kernel Kernel regression / RFF approximation [implemented]"+ , " spline B-spline / natural cubic regression [implemented]"+ , " quantile Quantile regression (τ-quantile, MM-IRLS) [implemented]"+ , " gam Generalized Additive Model (additive B-splines + Ridge) [implemented]"+ , " rf Random Forest regression (CART + bagging + feature subset) [implemented]"+ , " multireg Multi-output regression (wide CSV; linear/kernel-rbf) [implemented]"+ , " doe Orthogonal arrays (L_n) for experimental designs [implemented]"+ , " taguchi Taguchi method (SN ratio + factor effects + inner/outer) [implemented]"+ , ""+ , " help Show this message"+ , " --help, -h, help Same as 'help'"+ , ""+ , "Run 'hanalyze regress' (or invoke without a subcommand) to see regression-specific options."+ ]++futureSubcommands :: [(String, String)]+futureSubcommands =+ [+ ]++isFutureSubcommand :: String -> Bool+isFutureSubcommand c = c `elem` map fst futureSubcommands++stubMessage :: String -> String+stubMessage c = case lookup c futureSubcommands of+ Just msg -> msg+ Nothing -> "subcommand '" ++ c ++ "' is not yet implemented"++main :: IO ()+main = getArgs >>= dispatch++dispatch :: [String] -> IO ()+dispatch [] = putStrLn helpMsg+dispatch ("--help":_) = putStrLn helpMsg+dispatch ("-h":_) = putStrLn helpMsg+dispatch ("help":_) = putStrLn helpMsg+dispatch ("info":rest) = runInfoCmd rest+dispatch ("hist":rest) = runHistCmd rest+dispatch ("regress":rest) = runRegressCmd rest+dispatch ("doe":rest) = runDoeCmd rest+dispatch ("taguchi":rest) = runTaguchiCmd rest+dispatch ("ridge":rest) = runRidgeCmd rest+dispatch ("kernel":rest) = runKernelCmd rest+dispatch ("spline":rest) = runSplineCmd rest+dispatch ("quantile":rest) = runQuantileCmd rest+dispatch ("gam":rest) = runGAMCmd rest+dispatch ("rf":rest) = runRFCmd rest+dispatch ("clean":rest) = runCleanCmd rest+dispatch ("melt":rest) = runMeltCmd rest+dispatch ("regrid":rest) = runRegridCmd rest+dispatch ("multireg":rest) = runMultiRegCmd rest+dispatch (cmd:_) | isFutureSubcommand cmd = do+ hPutStrLn stderr $ "hanalyze: " ++ stubMessage cmd+ hPutStrLn stderr " Run 'hanalyze help' to see implemented subcommands."+dispatch args = runRegressCmd args -- legacy / bare++runRegressCmd :: [String] -> IO ()+runRegressCmd args = case parseArgs args of+ Left err -> hPutStrLn stderr err+ Right cfg -> runConfig cfg++-- ---------------------------------------------------------------------------+-- info subcommand+-- ---------------------------------------------------------------------------++runInfoCmd :: [String] -> IO ()+runInfoCmd args = do+ let (lopts, rest) = parseLoadOpts args+ case rest of+ [] -> hPutStrLn stderr+ "Usage: hanalyze info <file> [--no-header] [--skip N] [--comment CH] [--strict]"+ (file:_) -> do+ result <- loadAutoSafeWith lopts file+ case result of+ Left err -> hPutStrLn stderr ("Parse error: " ++ err)+ Right (df, lg) -> do+ Log.printLogReport lg+ printDataFrameInfo file df++-- | 共通フラグを切り出す: '--no-header' / '--skip N' / '--comment CH' /+-- '--strict' を 'LoadOpts' に集約し、残った位置引数を返す。+parseLoadOpts :: [String] -> (LoadOpts, [String])+parseLoadOpts = go defaultLoadOpts []+ where+ go acc rs [] = (acc, reverse rs)+ go acc rs ("--no-header":xs) = go acc { loNoHeader = True } rs xs+ go acc rs ("--strict":xs) = go acc { loStrict = True } rs xs+ go acc rs ("--no-sniff":xs) = go acc { loSniff = False } rs xs+ go acc rs ("--skip":n:xs)+ | Just k <- readMaybeInt n = go acc { loSkip = k } rs xs+ go acc rs ("--comment":cs:xs)+ | (c:_) <- cs = go acc { loComment = Just c } rs xs+ go acc rs (x:xs) = go acc (x:rs) xs++readMaybeInt :: String -> Maybe Int+readMaybeInt s = case reads s of+ [(n, "")] -> Just n+ _ -> Nothing++-- ---------------------------------------------------------------------------+-- 簡易プロファイリング (Phase 2)+-- ---------------------------------------------------------------------------++-- | アクションの実行時間を計測し、結果と合わせて (経過秒, 結果) を返す。+timed :: IO a -> IO (Double, a)+timed act = do+ t0 <- getCurrentTime+ r <- act+ t1 <- getCurrentTime+ return (realToFrac (diffUTCTime t1 t0), r)++-- | 1 段階のタイマー出力。" [Standardize] 12.3 ms" / " [Auto-HP] 58.21 s" 形式。+printPhase :: String -> Double -> IO ()+printPhase label sec+ | sec >= 1.0 = printf " [%-15s] %7.2f s\n" label sec+ | otherwise = printf " [%-15s] %7.0f ms\n" label (sec * 1000)++-- 未使用警告抑止+_dummyTime :: UTCTime -> UTCTime+_dummyTime = id++-- ---------------------------------------------------------------------------+-- clean subcommand (Phase C)+-- ---------------------------------------------------------------------------++runCleanCmd :: [String] -> IO ()+runCleanCmd args0 = do+ let (lopts, args1) = parseLoadOpts args0+ (rules, out, args2) = parseCleanFlags args1+ case args2 of+ [] -> hPutStrLn stderr cleanUsage+ (file:_) -> do+ result <- loadAutoSafeWith lopts file+ case result of+ Left err -> hPutStrLn stderr ("Parse error: " ++ err)+ Right (df0, lg0) -> do+ Log.printLogReport lg0+ let (df1, lg1) = Clean.cleanPipeline rules df0+ Log.printLogReport lg1+ case out of+ Nothing -> do+ -- 出力ファイル指定なし: info を出して終わり+ putStrLn "Cleaned DataFrame:"+ putStrLn $ " Rows / Cols: "+ ++ show (fst (DX.dimensions df1)) ++ " × "+ ++ show (length (DX.columnNames df1))+ putStrLn " Columns:"+ mapM_ (TIO.putStrLn . (" - " <>)) (DX.columnNames df1)+ Just path -> do+ -- TODO: 簡易 CSV 書出し。今は警告だけ出す。+ hPutStrLn stderr+ ("(--output " ++ path ++ " は未実装。ライブラリ API "+ ++ "Clean.cleanPipeline + Hackage writeCsv を直接お使いください)")++cleanUsage :: String+cleanUsage = unlines+ [ "Usage: hanalyze clean <file> [--rule COL=RULE]... [--output FILE] [load opts]"+ , ""+ , "Rules (各列に適用):"+ , " StripUnits \"12.3kg\" → 12.3"+ , " ParseCurrency \"$1,234.56\" → 1234.56"+ , " ParseDecimalEU \"3,14\" → 3.14 (decimal point が ,)"+ , " TrimText 前後空白を除去"+ , " CoerceNumeric 上記 3 種を順に試す万能変換"+ , ""+ , "例:"+ , " hanalyze clean data/raw.csv \\"+ , " --rule price=ParseCurrency \\"+ , " --rule weight=StripUnits \\"+ , " --rule note=TrimText"+ , ""+ , "Load opts: --no-header / --skip N / --comment CH / --delim CH / --strict / --no-sniff"+ ]++parseCleanFlags+ :: [String] -> ([(T.Text, Clean.ColumnRule)], Maybe FilePath, [String])+parseCleanFlags = go [] Nothing []+ where+ go rs out kept [] = (reverse rs, out, reverse kept)+ go rs out kept ("--rule":spec:xs) = case parseRuleSpec spec of+ Just r -> go (r:rs) out kept xs+ Nothing -> go rs out kept xs+ go rs _ kept ("--output":p:xs) = go rs (Just p) kept xs+ go rs _ kept ("-o":p:xs) = go rs (Just p) kept xs+ go rs out kept (x:xs) = go rs out (x:kept) xs++-- ---------------------------------------------------------------------------+-- melt subcommand (Phase B/C — wide → long)+-- ---------------------------------------------------------------------------++runMeltCmd :: [String] -> IO ()+runMeltCmd args0 = do+ let (lopts, args1) = parseLoadOpts args0+ (mopts, args2) = parseMeltFlags args1+ case args2 of+ [] -> hPutStrLn stderr meltUsage+ (file:_) -> case (moIds mopts, moVars mopts) of+ ([], _) -> hPutStrLn stderr "melt: --id COL1,COL2,... 必須"+ (_, []) -> hPutStrLn stderr "melt: --vars COL1,COL2,... 必須"+ (ids, vars) -> do+ result <- loadAutoSafeWith lopts file+ case result of+ Left err -> hPutStrLn stderr ("Parse error: " ++ err)+ Right (df0, lg) -> do+ Log.printLogReport lg+ let df1 = Pp.meltLonger ids vars+ (moVarName mopts) (moValueName mopts)+ (moParseVar mopts) df0+ (nrows, ncols) = DX.dimensions df1+ putStrLn "Long-form DataFrame:"+ putStrLn $ " Rows / Cols: " ++ show nrows ++ " × " ++ show ncols+ putStrLn " Columns:"+ mapM_ (TIO.putStrLn . (" - " <>)) (DX.columnNames df1)+ case moOut mopts of+ Just path -> do+ writeMeltedCsv path df1+ putStrLn $ "Wrote " ++ path+ Nothing -> return ()++-- | melt 結果を簡易 CSV (Hackage の writeCsv 経由) で書き出す。+writeMeltedCsv :: FilePath -> DXD.DataFrame -> IO ()+writeMeltedCsv path df = DX.writeCsv path df++data MeltOpts = MeltOpts+ { moIds :: [T.Text]+ , moVars :: [T.Text]+ , moVarName :: T.Text+ , moValueName :: T.Text+ , moParseVar :: Bool+ , moOut :: Maybe FilePath+ } deriving (Show)++defaultMeltOpts :: MeltOpts+defaultMeltOpts = MeltOpts [] [] "variable" "value" True Nothing++parseMeltFlags :: [String] -> (MeltOpts, [String])+parseMeltFlags = go defaultMeltOpts []+ where+ splitCSV = map T.pack . filter (not . null) . wordsBy (== ',')+ go acc kept [] = (acc, reverse kept)+ go acc kept ("--id":v:xs) = go acc { moIds = splitCSV v } kept xs+ go acc kept ("--vars":v:xs) = go acc { moVars = splitCSV v } kept xs+ go acc kept ("--var":v:xs) = go acc { moVarName = T.pack v } kept xs+ go acc kept ("--value":v:xs) = go acc { moValueName = T.pack v } kept xs+ go acc kept ("--no-parse-var":xs) = go acc { moParseVar = False } kept xs+ go acc kept ("--output":p:xs) = go acc { moOut = Just p } kept xs+ go acc kept ("-o":p:xs) = go acc { moOut = Just p } kept xs+ go acc kept (x:xs) = go acc (x:kept) xs++wordsBy :: (Char -> Bool) -> String -> [String]+wordsBy p s = case dropWhile p s of+ "" -> []+ s' -> let (w, rest) = break p s' in w : wordsBy p rest++-- ---------------------------------------------------------------------------+-- regrid subcommand (Phase G5)+-- ---------------------------------------------------------------------------++data RegridCliOpts = RegridCliOpts+ { rcId :: T.Text+ , rcZ :: T.Text+ , rcY :: T.Text+ , rcN :: Int+ , rcInterp :: Interp.InterpKind+ , rcGrid :: AG.GridKind+ , rcZBounds :: Pp.ZBoundsMode+ , rcReport :: Maybe FilePath+ , rcReportExtra :: Bool+ , rcOut :: Maybe FilePath+ } deriving (Show)++defaultRegridCliOpts :: RegridCliOpts+defaultRegridCliOpts = RegridCliOpts+ { rcId = "id"+ , rcZ = "z"+ , rcY = "y"+ , rcN = 30+ , rcInterp = Interp.PCHIP+ , rcGrid = AG.Adaptive+ , rcZBounds = Pp.ZIntersection+ , rcReport = Nothing+ , rcReportExtra = False+ , rcOut = Nothing+ }++parseRegridFlags :: [String] -> (RegridCliOpts, [String])+parseRegridFlags = go defaultRegridCliOpts []+ where+ go acc kept [] = (acc, reverse kept)+ go acc kept ("--id":v:xs) = go acc { rcId = T.pack v } kept xs+ go acc kept ("--z":v:xs) = go acc { rcZ = T.pack v } kept xs+ go acc kept ("--y":v:xs) = go acc { rcY = T.pack v } kept xs+ go acc kept ("--n":v:xs) =+ go acc { rcN = maybe 30 id (readMaybe v) } kept xs+ go acc kept ("--interp":v:xs) =+ let k = case v of+ "linear" -> Interp.Linear+ "spline" -> Interp.NaturalSpline+ "natural" -> Interp.NaturalSpline+ "naturalspline" -> Interp.NaturalSpline+ "pchip" -> Interp.PCHIP+ _ -> Interp.PCHIP+ in go acc { rcInterp = k } kept xs+ go acc kept ("--grid":v:xs) =+ let g = case v of "uniform" -> AG.Uniform+ "adaptive" -> AG.Adaptive+ _ -> AG.Adaptive+ in go acc { rcGrid = g } kept xs+ go acc kept ("--zrange":v:xs) =+ let z = case v of "intersect" -> Pp.ZIntersection+ "intersection" -> Pp.ZIntersection+ "union" -> Pp.ZUnion+ _ -> Pp.ZIntersection+ in go acc { rcZBounds = z } kept xs+ go acc kept ("--report":p:xs) = go acc { rcReport = Just p } kept xs+ go acc kept ("--report-extra":xs) = go acc { rcReportExtra = True } kept xs+ go acc kept ("--output":p:xs) = go acc { rcOut = Just p } kept xs+ go acc kept ("-o":p:xs) = go acc { rcOut = Just p } kept xs+ go acc kept (x:xs) = go acc (x:kept) xs++runRegridCmd :: [String] -> IO ()+runRegridCmd args0 = do+ let (lopts, args1) = parseLoadOpts args0+ (rcOpts, args2) = parseRegridFlags args1+ case args2 of+ [] -> hPutStrLn stderr regridUsage+ (file:_) -> do+ result <- loadAutoSafeWith lopts file+ case result of+ Left err -> hPutStrLn stderr ("Parse error: " ++ err)+ Right (df0, lg) -> do+ Log.printLogReport lg+ let opts = Pp.RegridOpts+ { Pp.roInterp = rcInterp rcOpts+ , Pp.roGridKind = rcGrid rcOpts+ , Pp.roN = rcN rcOpts+ , Pp.roZBoundsMode = rcZBounds rcOpts+ , Pp.roCoarseN = 200+ , Pp.roEpsRatio = 0.05+ }+ rr = Pp.regridLong (rcId rcOpts) (rcZ rcOpts) (rcY rcOpts)+ opts df0+ df1 = Pp.rrDataFrame rr+ (nrows, ncols) = DX.dimensions df1+ putStrLn "Regridded long-form DataFrame:"+ putStrLn $ " Rows / Cols: " ++ show nrows ++ " × " ++ show ncols+ putStrLn $ " Z range: ["+ ++ show (Pp.rrZMin rr) ++ ", "+ ++ show (Pp.rrZMax rr) ++ "]"+ putStrLn $ " N grid: " ++ show (length (Pp.rrZGrid rr))+ putStrLn $ " IDs: " ++ show (length (Pp.rrIds rr))+ case rcOut rcOpts of+ Just path -> do+ DX.writeCsv path df1+ putStrLn $ "Wrote " ++ path+ Nothing -> return ()+ case rcReport rcOpts of+ Just path -> do+ let kindStr = case rcInterp rcOpts of+ Interp.Linear -> "Linear"+ Interp.NaturalSpline -> "NaturalSpline"+ Interp.PCHIP -> "PCHIP"+ gridStr = case rcGrid rcOpts of+ AG.Uniform -> "Uniform"+ AG.Adaptive -> "Adaptive"+ zbStr = case rcZBounds rcOpts of+ Pp.ZIntersection -> "intersect"+ Pp.ZUnion -> "union"+ perObs = [ (i, pts) | (i, pts, _) <- Pp.rrPerIdInterp rr ]+ perInterp = [ (i, zip (Pp.rrZGrid rr) (map f (Pp.rrZGrid rr)))+ | (i, _, f) <- Pp.rrPerIdInterp rr ]+ perSummary = [ ( Pp.piId s, Pp.piNObserved s+ , Pp.piZMin s, Pp.piZMax s+ , Pp.piExtrapBelow s, Pp.piExtrapAbove s+ , Pp.piResidualMax s)+ | s <- Pp.rrPerIdStats rr ]+ perYRange = [ let ysOrig = map snd pts+ ysGrid = map (\(_,y) -> y) gys+ in (i, minimum ysOrig, maximum ysOrig+ , minimum ysGrid, maximum ysGrid)+ | ((i, pts, _), gys) <-+ zip (Pp.rrPerIdInterp rr) (map snd perInterp)+ ]+ ir = RB.InterpReport+ { RB.irTitle = "Regrid summary"+ , RB.irInterpKind = kindStr+ , RB.irGridKind = gridStr+ , RB.irN = rcN rcOpts+ , RB.irZBoundsMode = zbStr+ , RB.irZMin = Pp.rrZMin rr+ , RB.irZMax = Pp.rrZMax rr+ , RB.irPerIdObserved = perObs+ , RB.irPerIdInterpY = perInterp+ , RB.irGrid = Pp.rrZGrid rr+ , RB.irDensity = Pp.rrDensity rr+ , RB.irPerIdSummary = perSummary+ , RB.irExtraEnabled = rcReportExtra rcOpts+ , RB.irPerIdYRange = perYRange+ }+ RB.renderReport path+ (RB.defaultReportConfig "Regrid report")+ [RB.secInterpolation ir]+ putStrLn $ "Wrote report " ++ path+ Nothing -> return ()++regridUsage :: String+regridUsage = unlines+ [ "Usage: hanalyze regrid <file> [options] [load opts]"+ , ""+ , "歯抜けの long-form データ [id, z, y] を共通 grid に揃える。"+ , ""+ , " --id COL id 列名 (default: id)"+ , " --z COL z 列名 (default: z)"+ , " --y COL y 列名 (default: y)"+ , " --n N grid 点数 (default: 30)"+ , " --interp KIND linear | spline | pchip (default: pchip)"+ , " --grid KIND uniform | adaptive (default: adaptive)"+ , " --zrange MODE intersect | union (default: intersect)"+ , " --output FILE 揃った long-form を CSV で出力"+ , " --report FILE HTML レポート (R1-R7 必須要素)"+ , " --report-extra --report に R8-R10 オプション要素も追加"+ , ""+ , "例:"+ , " hanalyze regrid data/io/potential_long_jagged.csv \\"+ , " --id name --z z --y y --n 30 \\"+ , " --interp pchip --grid adaptive --zrange intersect \\"+ , " --output regridded.csv --report regrid.html --report-extra"+ ]++meltUsage :: String+meltUsage = unlines+ [ "Usage: hanalyze melt <file> --id COL1,COL2,... --vars COL1,COL2,..."+ , " [--var NAME] [--value NAME]"+ , " [--no-parse-var] [--output FILE] [load opts]"+ , ""+ , "wide-form CSV を long-form (tidy) に展開する。"+ , ""+ , " --id そのまま残す列 (例: name,x1,x2)"+ , " --vars 縦方向に展開する wide 列 (例: 1,2,3,4,5,6,7,8,9,10)"+ , " --var 新しい variable 列名 (default: 'variable'; 例: --var t)"+ , " --value 新しい value 列名 (default: 'value'; 例: --value y)"+ , " --no-parse-var variable 列を Double に parse せず Text のまま残す"+ , " --output FILE 結果を CSV として書き出す"+ , ""+ , "例:"+ , " hanalyze melt data/io/wide_sample.csv \\"+ , " --id name,x1,x2 \\"+ , " --vars 1,2,3,4,5,6,7,8,9,10 \\"+ , " --var t --value y \\"+ , " --output data/io/melted_sample.csv"+ ]++parseRuleSpec :: String -> Maybe (T.Text, Clean.ColumnRule)+parseRuleSpec s = case break (== '=') s of+ (col, '=':rule) | not (null col), not (null rule) ->+ case rule of+ "StripUnits" -> Just (T.pack col, Clean.StripUnits)+ "ParseCurrency" -> Just (T.pack col, Clean.ParseCurrency)+ "ParseDecimalEU" -> Just (T.pack col, Clean.ParseDecimalEU)+ "TrimText" -> Just (T.pack col, Clean.TrimText)+ "CoerceNumeric" -> Just (T.pack col, Clean.CoerceNumeric)+ _ -> Nothing+ _ -> Nothing++printDataFrameInfo :: FilePath -> DXD.DataFrame -> IO ()+printDataFrameInfo file df = do+ let n = (fst (DX.dimensions df))+ cols = DX.columnNames df+ putStrLn $ "File: " ++ file+ putStrLn $ "Rows: " ++ show n+ putStrLn $ "Columns: " ++ show (length cols)+ putStrLn ""+ printf " %-20s %-7s %5s %10s %10s %10s %10s %10s\n"+ ("name" :: String) ("type" :: String) ("n" :: String)+ ("min" :: String) ("max" :: String) ("mean" :: String)+ ("median" :: String) ("sd" :: String)+ putStrLn (replicate 92 '-')+ mapM_ (printColInfo df) cols++printColInfo :: DXD.DataFrame -> T.Text -> IO ()+printColInfo df name = case getDoubleVec name df of+ Just v -> do+ let xs = V.toList v+ m = length xs+ mn = if null xs then 0 else minimum xs+ mx = if null xs then 0 else maximum xs+ mean = if null xs then 0 else sum xs / fromIntegral m+ ss = sort xs+ med = if m == 0 then 0 else ss !! (m `div` 2)+ var = if m <= 1 then 0+ else sum [ (x - mean)^(2 :: Int) | x <- xs ]+ / fromIntegral (m - 1)+ sd_ = sqrt var+ printf " %-20s %-7s %5d %10.4f %10.4f %10.4f %10.4f %10.4f\n"+ (T.unpack name) ("numeric" :: String) m mn mx mean med sd_+ Nothing -> case getMaybeTextVec name df of+ Just v -> do+ let raw = V.toList v+ m = length raw+ xsOnly = [ x | Just x <- raw ]+ nMissNull = length [ () | Nothing <- raw ]+ nMissNA = length (filter Pp.isNAString xsOnly)+ nMiss = nMissNull + nMissNA+ uniq = Set.size (Set.fromList xsOnly)+ topN = take 3 (countTop xsOnly)+ topStr = intercalate ", "+ [ T.unpack k ++ "(" ++ show c ++ ")" | (k, c) <- topN ]+ missStr = if nMiss > 0+ then " NA=" ++ show nMiss+ else ""+ printf " %-20s %-7s %5d unique=%-3d top: %s%s\n"+ (T.unpack name) ("text" :: String) m uniq topStr missStr+ Nothing -> printf " %-20s %-7s ? (列の取り出しに失敗)\n"+ (T.unpack name) ("?" :: String)++-- | Count occurrences and return descending list.+countTop :: Ord a => [a] -> [(a, Int)]+countTop xs =+ let counts = foldr (\x -> insertWithInc x) [] xs+ insertWithInc x [] = [(x, 1)]+ insertWithInc x ((y, c) : rest)+ | x == y = (y, c + 1) : rest+ | otherwise = (y, c) : insertWithInc x rest+ sorted = qSortBy (\(_, a) (_, b) -> compare b a) counts+ in sorted+ where+ qSortBy _ [] = []+ qSortBy f (p:rs) = qSortBy f [x | x <- rs, f x p == LT || f x p == EQ]+ ++ [p]+ ++ qSortBy f [x | x <- rs, f x p == GT]++-- ---------------------------------------------------------------------------+-- hist subcommand+-- ---------------------------------------------------------------------------++runHistCmd :: [String] -> IO ()+runHistCmd args0 =+ let (lopts, args) = parseLoadOpts args0+ in case parseHistArgs args of+ Left err -> hPutStrLn stderr err+ Right ho -> runHistOpts (ho { hoLoadOpts = lopts })++data HistOpts = HistOpts+ { hoFile :: FilePath+ , hoCol :: T.Text+ , hoFit :: Maybe Distribution+ , hoFormat :: OutputFormat+ , hoOut :: FilePath+ , hoLoadOpts :: LoadOpts+ } deriving (Show)++parseHistArgs :: [String] -> Either String HistOpts+parseHistArgs (file : col : rest) = goHistOpts rest+ HistOpts { hoFile = file, hoCol = T.pack col+ , hoFit = Nothing, hoFormat = HTML, hoOut = "histogram.html"+ , hoLoadOpts = defaultLoadOpts }+parseHistArgs _ = Left $ unlines+ [ "Usage: hanalyze hist <file> <col> [options]"+ , ""+ , "Options:"+ , " --fit DIST [PARAMS...] overlay theoretical density"+ , " (e.g. --fit normal 0 1, --fit poisson 3)"+ , " --format html|png|svg output format (default: html)"+ , " --out FILE output file path (default: histogram.html)"+ ]++goHistOpts :: [String] -> HistOpts -> Either String HistOpts+goHistOpts [] ho = Right ho+goHistOpts ("--fit" : rest) ho = case rest of+ (name : rest') ->+ let (paramStrs, rest'') = span isNumericToken rest'+ params = map read paramStrs :: [Double]+ in case parseDistribution name params of+ Left err -> Left ("--fit: " ++ err)+ Right d -> goHistOpts rest'' (ho { hoFit = Just d })+ [] -> Left "--fit requires a distribution name (e.g. --fit normal 0 1)"+goHistOpts ("--format" : v : rest) ho =+ case parseFormat v of+ Left err -> Left err+ Right f -> goHistOpts rest (ho { hoFormat = f })+goHistOpts ("-f" : v : rest) ho = goHistOpts ("--format" : v : rest) ho+goHistOpts ("--out" : v : rest) ho = goHistOpts rest (ho { hoOut = v })+goHistOpts (flag : _) _ =+ Left ("hist: unexpected argument '" ++ flag ++ "' (try 'hanalyze hist' for usage)")++runHistOpts :: HistOpts -> IO ()+runHistOpts ho = do+ result <- loadAutoSafeWith (hoLoadOpts ho) (hoFile ho)+ case result of+ Left err -> hPutStrLn stderr ("Parse error: " ++ err)+ Right (df, lg) -> do+ Log.printLogReport lg+ case getDoubleVec (hoCol ho) df of+ Nothing -> hPutStrLn stderr $+ "Error: column '" ++ T.unpack (hoCol ho) ++ "' not found or not numeric"+ Just xVec -> do+ let vals = V.toList xVec+ histCfg = defaultConfig ("Histogram: " <> hoCol ho)+ outPath = hoOut ho+ fmt = hoFormat ho+ case hoFit ho of+ Nothing -> do+ histogramPlotFile fmt outPath histCfg (hoCol ho) vals Nothing+ putStrLn $ "Histogram: " ++ outPath+ Just dist -> do+ histogramWithDensityFile fmt outPath histCfg (hoCol ho) vals Nothing dist+ putStrLn $ "Histogram + density: " ++ outPath+ openInBrowser outPath++runConfig :: Config -> IO ()+runConfig cfg = do+ -- Warn: PI with non-Gaussian falls back to CI+ case (cfgBand cfg, cfgDist cfg, cfgModel cfg) of+ (PI _, fam, GLM) | fam /= Gaussian ->+ hPutStrLn stderr "Warning: PI is only exact for Gaussian. Using CI with same level."+ _ -> return ()+ -- Warn: GP only supports single x/y+ case cfgModel cfg of+ GP | length (cfgXCols cfg) /= 1 || length (cfgYCols cfg) /= 1 ->+ hPutStrLn stderr "Warning: GP requires exactly one x column and one y column."+ _ -> return ()++ result <- loadAutoSafeWith (cfgLoadOpts cfg) (cfgFile cfg)+ case result of+ Left err -> putStrLn ("Parse error: " ++ err)+ Right (df, lg) -> do+ Log.printLogReport lg+ putStrLn $ "Loaded " ++ show ((fst (DX.dimensions df))) ++ " rows from " ++ cfgFile cfg+ putStrLn "Columns:"+ mapM_ (TIO.putStrLn . (" - " <>)) (DX.columnNames df)++ let fmt = cfgFormat cfg+ xCol1 = head (cfgXCols cfg)++ -- ── Histogram mode ────────────────────────────────────────────────────+ if cfgHistMode cfg+ then runHistogram cfg df fmt xCol1+ else case cfgModel cfg of+ GP -> runGP cfg df xCol1+ HBM -> runHBM cfg df xCol1+ _ -> runAnalysis cfg df fmt xCol1++-- ---------------------------------------------------------------------------+-- Mixed model (LME / GLMM)+-- ---------------------------------------------------------------------------++runMixedModel :: Config -> DXD.DataFrame -> OutputFormat -> T.Text -> T.Text -> T.Text -> IO ()+runMixedModel cfg df fmt xCol1 yCol grpCol = do+ let colDegs = applyDegreeSpec (cfgDegree cfg) (cfgXCols cfg)+ (dist, lnk) = case cfgModel cfg of+ LM -> (Gaussian, Identity)+ GLM -> (cfgDist cfg, cfgLink cfg)+ _ -> (Gaussian, Identity) -- unreachable (NoReg/GP)++ let mResult = case cfgModel cfg of+ LM -> fitLMEDataFrame colDegs grpCol yCol df+ GLM -> fitGLMMDataFrame dist lnk colDegs grpCol yCol df+ _ -> Nothing++ case mResult of+ Nothing -> putStrLn "\nError: column(s) not found or not numeric/text"+ Just gr -> do+ let modelKind = case cfgModel cfg of+ LM -> "LME (Gaussian, exact EM)"+ GLM -> "GLMM (" ++ modelLabel dist lnk ++ ", Laplace)"+ _ -> ""+ cs = coeffList (glmmFixed gr)++ putStrLn $ "\nModel: " ++ T.unpack yCol ++ " ~ "+ ++ modelFormula colDegs+ ++ " [" ++ modelKind ++ " | group: " ++ T.unpack grpCol ++ "]"++ putStrLn "Fixed effects:"+ mapM_ (\(lbl, v) -> printf " %-30s = %9.4f\n" lbl v)+ (zip (multiCoeffLabels colDegs) cs)++ putStrLn "Variance components:"+ printf " %-30s = %9.4f\n" ("σ²_u (" ++ T.unpack grpCol ++ ")") (glmmRandVar gr)+ case cfgModel cfg of+ LM -> printf " %-30s = %9.4f\n" ("σ² (residual)" :: String) (glmmResidVar gr)+ _ -> printf " %-30s (fixed by family)\n" ("σ² (residual)" :: String)+ printf " %-30s = %9.4f (%d%% between-group)\n"+ ("ICC" :: String) (glmmICC gr)+ (round (glmmICC gr * 100) :: Int)++ putStrLn $ "BLUPs (" ++ T.unpack grpCol ++ "):"+ mapM_ (\(g, u) -> printf " %-12s = %+9.4f\n" g u)+ (zip (map T.unpack (V.toList (glmmGroups gr))) (V.toList (glmmBLUPs gr)))++ let suffix = " [" <> T.pack modelKind <> " | group: " <> grpCol <> "]"++ -- Scatter with group-level fitted lines (single x only)+ case (length (cfgXCols cfg), getDoubleVec xCol1 df, getDoubleVec yCol df, getTextVec grpCol df) of+ (1, Just xVec, Just yVec, Just gVec) -> do+ let ptData = zip3 (V.toList gVec) (V.toList xVec) (V.toList yVec)+ lnData = computeGroupLines lnk cs colDegs+ (glmmGroups gr) (glmmBLUPs gr) xVec+ scatterPath = "scatter.html"+ scatterCfg = defaultConfig (xCol1 <> " vs " <> yCol <> suffix)+ scatterWithGroupsFile fmt scatterPath scatterCfg xCol1 yCol ptData lnData+ putStrLn $ "\nScatter plot: " ++ scatterPath+ openInBrowser scatterPath+ _ ->+ putStrLn "\n(Scatter plot skipped for multiple x columns)"++ -- Predicted vs Actual+ case getDoubleVec yCol df of+ Nothing -> return ()+ Just yVec -> do+ let pvsaPath = "pvsa.html"+ pvsaCfg = defaultConfig ("Predicted vs Actual" <> suffix)+ predictedVsActualFile fmt pvsaPath pvsaCfg (V.toList yVec) (fittedList (glmmFixed gr))+ putStrLn $ "Predicted vs Actual: " ++ pvsaPath+ openInBrowser pvsaPath++ -- ── HTML レポート生成 ──────────────────────────────────────────────────+ case cfgReport cfg of+ Nothing -> return ()+ Just path -> do+ -- WAIC/LOO 計算 (--waic 指定時、Gaussian/Identity の LME のみ)+ mModelSel <-+ if cfgWAIC cfg && dist == Gaussian && lnk == Identity+ then case (getDoubleVec yCol df, getTextVec grpCol df) of+ (Just yVec, Just gVec) -> do+ let xVecPairs = [ (xv, deg) | (xc, deg) <- colDegs+ , Just xv <- [getDoubleVec xc df] ]+ case xVecPairs of+ [] -> return Nothing+ _ -> do+ let dm = multiPolyDesignMatrix xVecPairs+ y = LA.fromList (V.toList yVec)+ groupLabels = V.toList (glmmGroups gr)+ blupsList = V.toList (glmmBLUPs gr)+ blupMap = zip groupLabels blupsList+ offsets = [ maybe 0 id (lookup g blupMap)+ | g <- V.toList gVec ]+ nSamples = 1000+ gen <- createSystemRandom+ llMat <- lmePosteriorLogLiks+ dm y offsets (glmmFixed gr) nSamples gen+ let w = waic llMat+ l = loo llMat+ printf " WAIC=%.2f LOO=%.2f p_WAIC=%.2f k̂>0.7: %d件 (条件付き)\n"+ (waicValue w) (looValue l) (waicPwaic w) (looKHatBad l)+ return (Just (w, l))+ _ -> return Nothing+ else return Nothing++ let rbCfg = RB.defaultReportConfig+ (T.pack modelKind <> ": " <> yCol <> " | " <> grpCol)+ scatterPlots =+ case (length (cfgXCols cfg), getDoubleVec xCol1 df, getDoubleVec yCol df, getTextVec grpCol df) of+ (1, Just xVec, Just yVec, Just gVec) ->+ let ptData = zip3 (V.toList gVec) (V.toList xVec) (V.toList yVec)+ lnData = computeGroupLines lnk cs colDegs (glmmGroups gr) (glmmBLUPs gr) xVec+ scCfg = defaultConfig (xCol1 <> " vs " <> yCol <> suffix)+ in [NamedPlot "vl-scatter" "グループ別散布図"+ (scatterWithGroups scCfg xCol1 yCol ptData lnData)]+ _ -> []+ pvsaPlots =+ case getDoubleVec yCol df of+ Just yVec ->+ let pvCfg = defaultConfig ("Predicted vs Actual" <> suffix)+ in [NamedPlot "vl-pvsa" "Predicted vs Actual"+ (predictedVsActual pvCfg (V.toList yVec) (fittedList (glmmFixed gr)))]+ Nothing -> []+ plots = scatterPlots ++ pvsaPlots+ sections = cliMixedSections cfg df dist lnk colDegs grpCol gr mModelSel plots+ RB.renderReport path rbCfg sections+ putStrLn $ "Report: " ++ path+ maybeExportReportPlots cfg path plots+ openInBrowser path++-- ---------------------------------------------------------------------------+-- GLM regression (no random effects)+-- ---------------------------------------------------------------------------++runRegression :: Config -> DXD.DataFrame -> OutputFormat -> T.Text -> T.Text -> IO ()+runRegression cfg df fmt xCol1 yCol = do+ let colDegs = applyDegreeSpec (cfgDegree cfg) (cfgXCols cfg)+ (dist, lnk) = case cfgModel cfg of+ LM -> (Gaussian, Identity)+ GLM -> (cfgDist cfg, cfgLink cfg)+ _ -> (Gaussian, Identity) -- unreachable (NoReg/GP)++ case fitGLMWithSmooth dist lnk colDegs (cfgBand cfg) 200 df yCol of+ Nothing -> putStrLn "\nError: column(s) not found or not numeric"+ Just (res, mSmooth) -> do+ let cs = coeffList res+ eq = equationLabel dist lnk colDegs cs++ putStrLn $ "\nModel: " ++ T.unpack yCol ++ " ~ "+ ++ modelFormula colDegs+ ++ " [" ++ modelLabel dist lnk ++ "]"+ mapM_ (\(lbl, v) -> printf " %-30s = %9.4f\n" lbl v)+ (zip (multiCoeffLabels colDegs) cs)+ printf " %-30s = %9.4f\n" (r2Label dist) (rSquared1 res)++ let bandLabel = case cfgBand cfg of+ NoBand -> ""+ CI level -> ", " ++ show (round (level*100) :: Int) ++ "% CI"+ PI level -> ", " ++ show (round (level*100) :: Int) ++ "% PI"+ titleSuffix = " [" <> T.pack (modelLabel dist lnk) <> T.pack bandLabel <> "]"++ case mSmooth of+ Just sf -> do+ let scatterPath = "scatter.html"+ scatterCfg = defaultConfig (xCol1 <> " vs " <> yCol <> titleSuffix)+ scatterWithSmoothFile fmt scatterPath scatterCfg eq df xCol1 yCol sf+ putStrLn $ "\nScatter plot: " ++ scatterPath+ openInBrowser scatterPath+ Nothing ->+ putStrLn "\n(Scatter plot skipped for multiple x columns)"++ case getDoubleVec yCol df of+ Nothing -> return ()+ Just yVec -> do+ let pvsaPath = "pvsa.html"+ pvsaCfg = defaultConfig ("Predicted vs Actual " <> titleSuffix)+ predictedVsActualFile fmt pvsaPath pvsaCfg (V.toList yVec) (fittedList res)+ putStrLn $ "Predicted vs Actual: " ++ pvsaPath+ openInBrowser pvsaPath++ -- ── HTML レポート生成 ──────────────────────────────────────────────────+ case cfgReport cfg of+ Nothing -> return ()+ Just path -> do+ let rbCfg = RB.defaultReportConfig+ (T.pack (modelLabel dist lnk)+ <> ": " <> yCol <> " ~ "+ <> T.pack (modelFormula colDegs))+ pvsaPlots = case getDoubleVec yCol df of+ Just yVec ->+ let pvCfg = defaultConfig ("Predicted vs Actual" <> titleSuffix)+ in [NamedPlot "vl-pvsa" "Predicted vs Actual"+ (predictedVsActual pvCfg (V.toList yVec) (fittedList res))]+ Nothing -> []++ -- ── WAIC/LOO-CV 計算 (--waic が指定された場合) ────────────────────+ mModelSelect <-+ if not (cfgWAIC cfg)+ then return Nothing+ else case getDoubleVec yCol df of+ Nothing -> return Nothing+ Just yVec -> do+ let xVecPairs = [ (xv, deg)+ | (xc, deg) <- colDegs+ , Just xv <- [getDoubleVec xc df] ]+ case xVecPairs of+ [] -> return Nothing+ _ -> do+ let dm = multiPolyDesignMatrix xVecPairs+ y = LA.fromList (V.toList yVec)+ nSamples = 1000 :: Int+ gen <- createSystemRandom+ llMat <- case dist of+ Gaussian -> lmPosteriorLogLiks dm y res nSamples gen+ _ -> do+ let (_, fisherInv) = fitGLMFull dist lnk dm y+ glmPosteriorLogLiks dist lnk dm y fisherInv res nSamples gen+ let w = waic llMat+ l = loo llMat+ printf " WAIC=%.2f LOO=%.2f p_WAIC=%.2f k̂>0.7: %d件\n"+ (waicValue w) (looValue l) (waicPwaic w) (looKHatBad l)+ return (Just (w, l))++ let sections = cliRegressSections cfg df dist lnk colDegs res mSmooth+ mModelSelect pvsaPlots+ RB.renderReport path rbCfg sections+ putStrLn $ "Report: " ++ path+ maybeExportReportPlots cfg path pvsaPlots+ openInBrowser path++-- ---------------------------------------------------------------------------+-- Regression / scatter dispatch (non-histogram path)+-- ---------------------------------------------------------------------------++runAnalysis :: Config -> DXD.DataFrame -> OutputFormat -> T.Text -> IO ()+runAnalysis cfg df fmt xCol1 = do+ let yCols = cfgYCols cfg+ effModel = if length yCols > 1 then NoReg+ else if cfgModel cfg == GP then NoReg -- GP はここに来ない+ else cfgModel cfg++ case effModel of+ -- ── No regression: scatter plot only ──────────────────────────────────+ NoReg ->+ case yCols of+ [yCol] -> do+ let scatterPath = "scatter.html"+ scatterCfg = defaultConfig (xCol1 <> " vs " <> yCol)+ scatterPlotFile fmt scatterPath scatterCfg df xCol1 yCol+ putStrLn $ "\nScatter plot: " ++ scatterPath+ openInBrowser scatterPath++ _ -> do+ let scatterPath = "scatter.html"+ scatterCfg = defaultConfig (xCol1 <> " vs " <> T.intercalate ", " yCols)+ scatterMultiYFile fmt scatterPath scatterCfg df xCol1 yCols+ putStrLn $ "\nScatter plot (multi-y): " ++ scatterPath+ openInBrowser scatterPath++ -- ── Regression (LM / GLM) ─────────────────────────────────────────────+ _ -> case yCols of+ [yCol] ->+ case cfgGroup cfg of+ Just grpCol -> runMixedModel cfg df fmt xCol1 yCol grpCol+ Nothing -> runRegression cfg df fmt xCol1 yCol+ _ -> do+ putStrLn "\nNote: regression with multiple y columns not supported. Plotting scatter only."+ let scatterPath = "scatter.html"+ scatterCfg = defaultConfig (xCol1 <> " vs " <> T.intercalate ", " yCols)+ scatterMultiYFile fmt scatterPath scatterCfg df xCol1 yCols+ putStrLn $ "Scatter plot (multi-y): " ++ scatterPath+ openInBrowser scatterPath++-- ---------------------------------------------------------------------------+-- GP regression+-- ---------------------------------------------------------------------------++runGP :: Config -> DXD.DataFrame -> T.Text -> IO ()+runGP cfg df xCol1 = do+ let yCol = head (cfgYCols cfg)+ case (getDoubleVec xCol1 df, getDoubleVec yCol df) of+ (Just xVec, Just yVec) -> do+ let xs = V.toList xVec+ ys = V.toList yVec+ p0 = initParamsFromData xs ys++ putStrLn "\nFitting GP kernels (this may take a moment)..."++ let kernelDefs = [(RBF, "RBF"), (Matern52, "Mat\xe9rn5/2"), (Periodic, "Periodic")]+ xMin = V.minimum xVec+ xMax = V.maximum xVec+ span' = max 1e-8 (xMax - xMin)+ testXs = [ xMin + fromIntegral i * span' / 199 | i <- [0 .. 199 :: Int] ]++ kfits <- mapM (\(ker, lbl) -> do+ putStrLn $ " Optimizing " ++ lbl ++ " ..."+ let params = optimizeGP ker xs ys p0+ model = GPModel ker params+ res = fitGP model xs ys testXs+ lml = logMarginalLikelihood xs ys ker params+ pd = gpPredData model xs ys+ return GPKernelFit+ { gkLabel = T.pack lbl+ , gkKernel = ker+ , gkParams = params+ , gkResult = res+ , gkLML = lml+ , gkPredData = pd+ }+ ) kernelDefs++ -- LML 降順にソート+ let sorted = foldr insertByLML [] kfits+ insertByLML x [] = [x]+ insertByLML x (y:ys') = if gkLML x >= gkLML y then x:y:ys'+ else y : insertByLML x ys'+ path = maybe "report.html" id (cfgReport cfg)+ rbCfg = RB.defaultReportConfig+ ("GP Regression: " <> xCol1 <> " \x2192 " <> yCol)+ sections = cliGPSections xCol1 yCol df xs ys testXs sorted++ RB.renderReport path rbCfg sections+ putStrLn $ "Report: " ++ path+ maybeExportReportPlots cfg path []+ openInBrowser path++ _ -> putStrLn "\nError: column(s) not found or not numeric"++-- ---------------------------------------------------------------------------+-- HBM (Bayesian linear regression via NUTS)+-- ---------------------------------------------------------------------------++runHBM :: Config -> DXD.DataFrame -> T.Text -> IO ()+runHBM cfg df xCol = do+ let yCols = cfgYCols cfg+ xCols = cfgXCols cfg+ case (yCols, xCols, getDoubleVec xCol df) of+ ([yCol], [_], Just xVec) ->+ case getDoubleVec yCol df of+ Nothing -> putStrLn $ "Error: y column '" ++ T.unpack yCol ++ "' not numeric"+ Just yVec -> do+ let xs = V.toList xVec+ ys = V.toList yVec+ putStrLn ""+ putStrLn "=== HBM Bayesian Linear Regression ==="+ printf " y = α + β·x + ε, α,β ~ Normal(0,10), ε ~ Normal(0,σ), σ ~ Exp(1)\n"+ printf " サンプリング: NUTS (AD 勾配 + dual averaging)\n"+ printf " N = %d 観測, x = %s, y = %s\n\n"+ (length xs) (T.unpack xCol) (T.unpack yCol)+ runHBMRegression xs ys xCol yCol df cfg+ _ ->+ putStrLn "Error: HBM requires exactly one x and one y column (numeric)"++runHBMRegression+ :: [Double] -> [Double] -> T.Text -> T.Text -> DXD.DataFrame -> Config -> IO ()+runHBMRegression xs ys xCol yCol df cfg = do+ let nutsCfg = HBMnuts.defaultNUTSConfig+ { HBMnuts.nutsIterations = 1500+ , HBMnuts.nutsBurnIn = 500+ , HBMnuts.nutsStepSize = 0.05+ }+ initP = Map.fromList+ [ ("alpha", 0.0), ("beta", 0.0), ("sigma", 1.0) ]+ hbmModel :: HBMod.ModelP ()+ hbmModel = do+ a <- HBMod.sample "alpha" (HBMod.Normal 0 10)+ b <- HBMod.sample "beta" (HBMod.Normal 0 10)+ s <- HBMod.sample "sigma" (HBMod.Exponential 1)+ mapM_ (\(x, y) ->+ let xC = realToFrac x+ in HBMod.observe "y" (HBMod.Normal (a + b * xC) s) [y])+ (zip xs ys)++ gen <- createSystemRandom+ chain <- HBMnuts.nuts hbmModel nutsCfg initP gen+ let acc = MCMCcore.acceptanceRate chain+ n = length (MCMCcore.chainSamples chain)+ printf " 受容率: %.1f%%, サンプル数: %d\n" (acc * 100 :: Double) n++ let aMean = maybe 0 id (MCMCcore.posteriorMean "alpha" chain)+ aSD = maybe 0 id (MCMCcore.posteriorSD "alpha" chain)+ bMean = maybe 0 id (MCMCcore.posteriorMean "beta" chain)+ bSD = maybe 0 id (MCMCcore.posteriorSD "beta" chain)+ sMean = maybe 0 id (MCMCcore.posteriorMean "sigma" chain)+ sSD = maybe 0 id (MCMCcore.posteriorSD "sigma" chain)+ printf " α = %+.4f ± %.4f\n" aMean aSD+ printf " β = %+.4f ± %.4f\n" bMean bSD+ printf " σ = %+.4f ± %.4f\n" sMean sSD++ case cfgReport cfg of+ Nothing -> return ()+ Just path -> do+ let smooth = makeSmooth xs chain+ fitted = [aMean + bMean * x | x <- xs]+ resid = zipWith (-) ys fitted+ yBar = sum ys / fromIntegral (length ys)+ tss = sum [(y - yBar) ^ (2::Int) | y <- ys]+ rss = sum [r ^ (2::Int) | r <- resid]+ r2 = if tss < 1e-12 then 0 else 1 - rss / tss++ mWaicLoo <-+ if cfgWAIC cfg+ then do+ let llMat = [ HBMod.perObsLogLiks hbmModel ps+ | ps <- MCMCcore.chainSamples chain ]+ w = waic llMat+ l = loo llMat+ printf " WAIC=%.2f LOO=%.2f p_WAIC=%.2f k̂>0.7: %d件\n"+ (waicValue w) (looValue l) (waicPwaic w) (looKHatBad l)+ return (Just (w, l))+ else return Nothing++ let mGraph = Just (Hanalyze.Viz.ModelGraph.buildMermaid (HBMod.buildModelGraph hbmModel))+ rbCfg = RB.defaultReportConfig+ ("HBM Linear Regression: " <> yCol <> " ~ " <> xCol)+ sections = cliHBMSections xCol yCol df xs ys chain mGraph mWaicLoo []+ RB.renderReport path rbCfg sections+ putStrLn $ "Report: " ++ path+ maybeExportReportPlots cfg path []+ openInBrowser path+ where+ -- 信用区間付き予測曲線: 各事後サンプルから μ* = α + β·x* を計算 → 分位点+ makeSmooth :: [Double] -> MCMCcore.Chain -> SmoothData+ makeSmooth xs0 ch =+ let alphas = MCMCcore.chainVals "alpha" ch+ betas = MCMCcore.chainVals "beta" ch+ xMin = minimum xs0+ xMax = maximum xs0+ ext = (xMax - xMin) * 0.5+ grid = [xMin - ext + i * (xMax - xMin + 2 * ext) / 99 | i <- [0..99]]+ atX x = let ss = sortListAsc (zipWith (\a b -> a + b * x) alphas betas)+ sn = length ss+ qAt p = ss !! min (sn-1) (max 0 (floor (p * fromIntegral sn) :: Int))+ in (qAt 0.5, qAt 0.025, qAt 0.975)+ (yMid, yLo, yHi) = unzip3 (map atX grid)+ in SmoothData+ { sdXs = grid, sdYs = yMid, sdLower = yLo, sdUpper = yHi+ , sdHasBand = True+ }++ sortListAsc :: [Double] -> [Double]+ sortListAsc = qs+ where+ qs [] = []+ qs (p:rs) = qs [x | x <- rs, x <= p] ++ [p] ++ qs [x | x <- rs, x > p]++-- ---------------------------------------------------------------------------+-- Histogram mode+-- ---------------------------------------------------------------------------++runHistogram :: Config -> DXD.DataFrame -> OutputFormat -> T.Text -> IO ()+runHistogram cfg df fmt xCol =+ case getDoubleVec xCol df of+ Nothing ->+ putStrLn $ "Error: column '" ++ T.unpack xCol ++ "' not found or not numeric"+ Just xVec -> do+ let vals = V.toList xVec+ histPath = "histogram.html"+ histCfg = defaultConfig ("Histogram: " <> xCol)+ case cfgFitDist cfg of+ Nothing -> do+ histogramPlotFile fmt histPath histCfg xCol vals Nothing+ putStrLn $ "\nHistogram: " ++ histPath+ Just dist -> do+ histogramWithDensityFile fmt histPath histCfg xCol vals Nothing dist+ putStrLn $ "\nHistogram + density: " ++ histPath+ openInBrowser histPath++-- ---------------------------------------------------------------------------+-- Group-level prediction helpers+-- ---------------------------------------------------------------------------++invLink :: LinkFn -> Double -> Double+invLink Identity eta = eta+invLink Log eta = exp eta+invLink Logit eta = 1.0 / (1.0 + exp (negate eta))+invLink Sqrt eta = eta * eta++-- | Generate per-group conditional fitted lines for visualization.+-- Only produces data when there is exactly one x column (scatter plot is 2D).+-- Returns [(group, xGrid, ŷ)] evaluated on a 100-point grid over [min(x), max(x)].+computeGroupLines+ :: LinkFn+ -> [Double] -- fixed coefficients [β₀, β₁, ..., βd]+ -> [(T.Text, Int)] -- x column / degree specs (length 1 → draw lines)+ -> V.Vector T.Text -- group labels (sorted)+ -> V.Vector Double -- BLUPs (same order as group labels)+ -> V.Vector Double -- observed x values (used to determine grid range)+ -> [(T.Text, Double, Double)]+computeGroupLines lnk coeffs colDegs groups blups xVec =+ case colDegs of+ [(_, deg)] | not (V.null xVec) ->+ let xMin = V.minimum xVec+ xMax = V.maximum xVec+ nGrid = 100 :: Int+ grid = [ xMin + fromIntegral i * (xMax - xMin) / fromIntegral (nGrid - 1)+ | i <- [0 .. nGrid - 1] ]+ b0 = head coeffs+ bs = tail coeffs+ etaAt x = b0 + sum (zipWith (*) bs [x ^ k | k <- [1 .. deg :: Int]])+ in [ (grp, x, invLink lnk (etaAt x + u))+ | (grp, u) <- zip (V.toList groups) (V.toList blups)+ , x <- grid ]+ _ -> []++-- ---------------------------------------------------------------------------+-- Formatting helpers+-- ---------------------------------------------------------------------------++modelLabel :: Family -> LinkFn -> String+modelLabel dist lnk = show dist ++ "/" ++ show lnk++r2Label :: Family -> String+r2Label Gaussian = "R²"+r2Label _ = "McFadden R²"++modelFormula :: [(T.Text, Int)] -> String+modelFormula colDegs = intercalate " + " (concatMap terms colDegs)+ where+ terms (col, deg) =+ [ T.unpack col ++ if k == 1 then "" else "^" ++ show k+ | k <- [1 .. deg]+ ]++multiCoeffLabels :: [(T.Text, Int)] -> [String]+multiCoeffLabels colDegs = "β₀ (intercept)" : zipWith fmt [1..] rest+ where+ rest = concatMap expand colDegs+ expand (col, deg) = [(col, k) | k <- [1 .. deg]]+ fmt i (col, k) =+ "β" ++ show (i :: Int) ++ " ("+ ++ T.unpack col+ ++ (if k == 1 then "" else "^" ++ show k)+ ++ ")"++-- | Generate a human-readable regression equation for single x-column models.+equationLabel :: Family -> LinkFn -> [(T.Text, Int)] -> [Double] -> Maybe T.Text+equationLabel _ _ colDegs _ | length colDegs /= 1 = Nothing+equationLabel _ _ _ coeffs | null coeffs = Nothing+equationLabel fam lnk [(col, deg)] coeffs = Just (T.pack label)+ where+ lhs = case (fam, lnk) of+ (Gaussian, Identity) -> "y"+ (_, Identity) -> "E[y]"+ _ -> show lnk ++ "(y)"++ b0 = head coeffs+ betas = tail coeffs++ termStr b k =+ let sign = if b >= 0 then " + " else " - "+ xStr = T.unpack col ++ if k == 1 then "" else "^" ++ show k+ in sign ++ printf "%.4f" (abs b :: Double) ++ xStr++ label = lhs ++ " = " ++ printf "%.4f" b0+ ++ concat (zipWith termStr betas [1 .. deg])+equationLabel _ _ _ _ = Nothing++-- ---------------------------------------------------------------------------+-- doe subcommand (Phase E1: orthogonal arrays)+-- ---------------------------------------------------------------------------++doeUsage :: String+doeUsage = unlines+ [ "Usage: hanalyze doe <action> [args...]"+ , ""+ , "Actions:"+ , " list List available standard arrays"+ , " ortho <NAME> [opts] Output an orthogonal array (L4/L8/L9/L12/L16/L18)"+ , ""+ , "ortho options:"+ , " -f, --factor NAME=v1,v2,... Assign a factor with comma-separated levels"+ , " (repeat for multiple factors; left-to-right = column 1, 2, ...)"+ , " --csv | --tsv | --pretty Output format (default: pretty)"+ , " --out FILE Write to file instead of stdout"+ , ""+ , "Examples:"+ , " hanalyze doe list"+ , " hanalyze doe ortho L9 --pretty"+ , " hanalyze doe ortho L9 -f temp=150,180,210 -f time=10,20,30 -f catalyst=A,B,C --csv"+ , " hanalyze doe ortho L8 -f A=low,high -f B=0,1 --out design.tsv --tsv"+ ]++runDoeCmd :: [String] -> IO ()+runDoeCmd [] = putStrLn doeUsage+runDoeCmd ["help"] = putStrLn doeUsage+runDoeCmd ["--help"] = putStrLn doeUsage+runDoeCmd ("list":_) = runDoeList+runDoeCmd ("ortho":rest) = runDoeOrtho rest+runDoeCmd (action:_) =+ hPutStrLn stderr ("doe: unknown action '" ++ action ++ "'\n" ++ doeUsage)++runDoeList :: IO ()+runDoeList = do+ putStrLn "Available standard orthogonal arrays:"+ mapM_ (\(name, descr) ->+ printf " %-16s %s\n" (T.unpack name) (T.unpack descr))+ OA.listArrays+ putStrLn ""+ putStrLn "Use 'hanalyze doe ortho <NAME>' to output a specific array."++data OrthoOpts = OrthoOpts+ { ooFactors :: [(T.Text, [T.Text])] -- name → comma-split levels+ , ooFormat :: OrthoOutFormat+ , ooOut :: Maybe FilePath+ } deriving (Show)++data OrthoOutFormat = OrthoCSV | OrthoTSV | OrthoPretty deriving (Show, Eq)++defaultOrthoOpts :: OrthoOpts+defaultOrthoOpts = OrthoOpts [] OrthoPretty Nothing++runDoeOrtho :: [String] -> IO ()+runDoeOrtho [] = hPutStrLn stderr ("doe ortho: missing array name\n" ++ doeUsage)+runDoeOrtho (nameStr : rest) =+ case OA.lookupOA (T.pack nameStr) of+ Nothing -> hPutStrLn stderr $+ "doe ortho: unknown array '" ++ nameStr+ ++ "' (try 'hanalyze doe list')"+ Just oa -> case parseOrthoOpts rest defaultOrthoOpts of+ Left err -> hPutStrLn stderr ("doe ortho: " ++ err)+ Right opts -> emitOrtho oa opts++parseOrthoOpts :: [String] -> OrthoOpts -> Either String OrthoOpts+parseOrthoOpts [] acc = Right acc+parseOrthoOpts (flag : rest) acc+ | flag `elem` ["-f", "--factor"] = case rest of+ (v : rest') -> case parseFactorSpec v of+ Left err -> Left err+ Right fac -> parseOrthoOpts rest' (acc { ooFactors = ooFactors acc ++ [fac] })+ [] -> Left "-f/--factor requires an argument like NAME=v1,v2,..."+ | flag == "--csv" = parseOrthoOpts rest (acc { ooFormat = OrthoCSV })+ | flag == "--tsv" = parseOrthoOpts rest (acc { ooFormat = OrthoTSV })+ | flag == "--pretty" = parseOrthoOpts rest (acc { ooFormat = OrthoPretty })+ | flag == "--out" = case rest of+ (v : rest') -> parseOrthoOpts rest' (acc { ooOut = Just v })+ [] -> Left "--out requires a file path"+ | otherwise = Left ("unexpected argument '" ++ flag ++ "'")++parseFactorSpec :: String -> Either String (T.Text, [T.Text])+parseFactorSpec s =+ case break (== '=') s of+ (name, '=' : levelsStr) | not (null name), not (null levelsStr) ->+ let levels = filter (not . T.null) (T.splitOn "," (T.pack levelsStr))+ in if null levels+ then Left ("factor '" ++ name ++ "' has no levels (use NAME=v1,v2,...)")+ else Right (T.pack name, levels)+ _ -> Left ("invalid factor spec '" ++ s ++ "' (expected NAME=v1,v2,...)")++emitOrtho :: OA.OA -> OrthoOpts -> IO ()+emitOrtho oa opts =+ case ooFactors opts of+ [] -> emitText (renderRaw (ooFormat opts) oa) (ooOut opts)+ fs -> do+ let specs = [ OA.FactorSpec name (map toLevelValue levels)+ | (name, levels) <- fs ]+ case OA.assignFactors oa specs of+ Left err -> hPutStrLn stderr ("doe ortho: " ++ T.unpack err)+ Right ad -> emitText (renderAssigned (ooFormat opts) ad) (ooOut opts)++toLevelValue :: T.Text -> OA.LevelValue+toLevelValue t = case reads (T.unpack t) :: [(Double, String)] of+ [(d, "")] -> OA.LNumeric d+ _ -> OA.LText t++renderRaw :: OrthoOutFormat -> OA.OA -> T.Text+renderRaw OrthoCSV = OA.renderRawCSV+renderRaw OrthoTSV = OA.renderRawTSV+renderRaw OrthoPretty = OA.renderRawPretty++renderAssigned :: OrthoOutFormat -> OA.AssignedDesign -> T.Text+renderAssigned OrthoCSV = OA.renderCSV+renderAssigned OrthoTSV = OA.renderTSV+renderAssigned OrthoPretty = OA.renderPretty++emitText :: T.Text -> Maybe FilePath -> IO ()+emitText txt Nothing = TIO.putStrLn txt+emitText txt (Just path) = do+ TIO.writeFile path txt+ putStrLn $ "Written: " ++ path++-- ---------------------------------------------------------------------------+-- taguchi subcommand (Phase E2: SN ratio + factor effects + inner/outer)+-- ---------------------------------------------------------------------------++taguchiUsage :: String+taguchiUsage = unlines+ [ "Usage: hanalyze taguchi <action> [args...]"+ , ""+ , "Actions:"+ , " sn <type> <values...> Compute a single SN ratio (dB)"+ , " type: smaller | larger | nominal | nominal-target=M"+ , ""+ , " analyze <ARRAY> -f F=v1,v2,... [-f ...] --csv FILE [--sntype TYPE] [--report [FILE]]"+ , " Analyze observations from a CSV file:"+ , " rows = inner runs, cols (after factor cols) = repetitions/outer."+ , " Computes per-row SN ratio, factor effects, and optimum levels."+ , " --report writes an interactive HTML report (default: taguchi.html)."+ , ""+ , " cross <INNER> <OUTER>"+ , " -f Fc=v1,v2,... [-f ...] Inner control factors"+ , " --noise Fn=v1,v2,... [...] Outer noise factors"+ , " [--out FILE] Output the cross-design CSV template"+ , ""+ , "SN types:"+ , " smaller smaller-the-better (e.g. defect rate)"+ , " larger larger-the-better (e.g. strength)"+ , " nominal nominal-the-best (mean^2 / variance)"+ , " nominal-target=M nominal with target value M"+ , ""+ , "Examples:"+ , " hanalyze taguchi sn smaller 1.2 1.5 0.9 1.1"+ , " hanalyze taguchi analyze L9 -f temp=150,180,210 -f time=10,20,30 -f cat=A,B,C"+ , " --csv runs.csv --sntype smaller"+ , " hanalyze taguchi cross L9 L4 -f temp=150,180,210 -f time=10,20,30 -f cat=A,B,C"+ , " --noise humidity=low,high --noise vibration=on,off --out cross.csv"+ ]++runTaguchiCmd :: [String] -> IO ()+runTaguchiCmd [] = putStrLn taguchiUsage+runTaguchiCmd ["help"] = putStrLn taguchiUsage+runTaguchiCmd ["--help"] = putStrLn taguchiUsage+runTaguchiCmd ("sn":rest) = runTaguchiSN rest+runTaguchiCmd ("analyze":rest) = runTaguchiAnalyze rest+runTaguchiCmd ("cross":rest) = runTaguchiCross rest+runTaguchiCmd (action:_) =+ hPutStrLn stderr ("taguchi: unknown action '" ++ action ++ "'\n" ++ taguchiUsage)++-- ── sn ──────────────────────────────────────────────────────────────────++runTaguchiSN :: [String] -> IO ()+runTaguchiSN [] = hPutStrLn stderr "taguchi sn: missing type and values"+runTaguchiSN (typeStr : valStrs)+ | null valStrs = hPutStrLn stderr "taguchi sn: need at least one value"+ | otherwise = case parseSNType typeStr of+ Left err -> hPutStrLn stderr ("taguchi sn: " ++ err)+ Right t ->+ let vals = mapM readMaybeD valStrs+ in case vals of+ Nothing -> hPutStrLn stderr "taguchi sn: non-numeric value(s)"+ Just xs -> do+ let eta = TG.snRatio t xs+ printf "SN(%s) = %.4f dB (n=%d)\n"+ (T.unpack (TG.snTypeName t)) eta (length xs)++parseSNType :: String -> Either String TG.SNType+parseSNType s = case s of+ "smaller" -> Right TG.SmallerBetter+ "smaller-better" -> Right TG.SmallerBetter+ "larger" -> Right TG.LargerBetter+ "larger-better" -> Right TG.LargerBetter+ "nominal" -> Right TG.NominalBest+ "nominal-best" -> Right TG.NominalBest+ _ | "nominal-target=" `isPrefixOfStr` s ->+ case readMaybeD (drop (length ("nominal-target=" :: String)) s) of+ Just m -> Right (TG.NominalBestTarget m)+ Nothing -> Left ("invalid target value in '" ++ s ++ "'")+ _ -> Left ("unknown SN type '" ++ s+ ++ "' (try smaller | larger | nominal | nominal-target=M)")++isPrefixOfStr :: String -> String -> Bool+isPrefixOfStr p s = take (length p) s == p++readMaybeD :: String -> Maybe Double+readMaybeD s = case reads s :: [(Double, String)] of+ [(v, "")] -> Just v+ _ -> Nothing++-- ── analyze ─────────────────────────────────────────────────────────────++data TgAnalyzeOpts = TgAnalyzeOpts+ { toFactors :: [(T.Text, [T.Text])]+ , toCSV :: Maybe FilePath+ , toSN :: TG.SNType+ , toReport :: Maybe FilePath+ } deriving (Show)++defaultTgAnalyzeOpts :: TgAnalyzeOpts+defaultTgAnalyzeOpts = TgAnalyzeOpts [] Nothing TG.SmallerBetter Nothing++runTaguchiAnalyze :: [String] -> IO ()+runTaguchiAnalyze args0 =+ let (lopts, args) = parseLoadOpts args0+ in case args of+ [] -> hPutStrLn stderr "taguchi analyze: missing array name"+ (arrayStr : rest) ->+ case OA.lookupOA (T.pack arrayStr) of+ Nothing -> hPutStrLn stderr $+ "taguchi analyze: unknown array '" ++ arrayStr ++ "'"+ Just oa -> case parseTgAnalyzeOpts rest defaultTgAnalyzeOpts of+ Left err -> hPutStrLn stderr ("taguchi analyze: " ++ err)+ Right opts -> case toCSV opts of+ Nothing -> hPutStrLn stderr "taguchi analyze: --csv FILE required"+ Just path -> doTaguchiAnalyze oa opts path lopts++parseTgAnalyzeOpts :: [String] -> TgAnalyzeOpts -> Either String TgAnalyzeOpts+parseTgAnalyzeOpts [] acc = Right acc+parseTgAnalyzeOpts (flag : rest) acc+ | flag `elem` ["-f", "--factor"] = case rest of+ (v : rs) -> case parseFactorSpec v of+ Left err -> Left err+ Right fac -> parseTgAnalyzeOpts rs+ (acc { toFactors = toFactors acc ++ [fac] })+ [] -> Left "-f/--factor requires NAME=v1,v2,..."+ | flag == "--csv" = case rest of+ (v : rs) -> parseTgAnalyzeOpts rs (acc { toCSV = Just v })+ [] -> Left "--csv requires a file path"+ | flag == "--sntype" = case rest of+ (v : rs) -> case parseSNType v of+ Left err -> Left err+ Right t -> parseTgAnalyzeOpts rs (acc { toSN = t })+ [] -> Left "--sntype requires an argument"+ | flag == "--report" = case rest of+ (v : rs) | not (null v) && head v /= '-' ->+ parseTgAnalyzeOpts rs (acc { toReport = Just v })+ _ -> parseTgAnalyzeOpts rest (acc { toReport = Just "taguchi.html" })+ | otherwise = Left ("unexpected argument '" ++ flag ++ "'")++doTaguchiAnalyze :: OA.OA -> TgAnalyzeOpts -> FilePath -> LoadOpts -> IO ()+doTaguchiAnalyze oa opts path lopts = do+ result <- loadAutoSafeWith lopts path+ case result of+ Left err -> hPutStrLn stderr ("Parse error: " ++ err)+ Right (df, lg) -> do+ Log.printLogReport lg+ let specs = [ OA.FactorSpec name (map toLevelValue lvls)+ | (name, lvls) <- toFactors opts ]+ case OA.assignFactors oa specs of+ Left err -> hPutStrLn stderr (T.unpack err)+ Right ad -> runAnalyzeWith ad opts df++runAnalyzeWith :: OA.AssignedDesign -> TgAnalyzeOpts -> DXD.DataFrame -> IO ()+runAnalyzeWith ad opts df = do+ let factorNames = map OA.fsName (OA.adFactors ad)+ yCols = filter (\c -> not (c `elem` factorNames) && c /= "Run")+ (DX.columnNames df)+ n = length (OA.adRows ad)+ when ((fst (DX.dimensions df)) /= n) $+ hPutStrLn stderr $+ "Warning: CSV has " ++ show ((fst (DX.dimensions df)))+ ++ " rows, expected " ++ show n+ if null yCols+ then hPutStrLn stderr+ "taguchi analyze: no observation columns found in CSV"+ else do+ -- Per-inner-run observations (skip non-numeric rows)+ let yMatrix =+ [ [ case getDoubleVec c df of+ Just v | i < V.length v -> v V.! i+ _ -> 0+ | c <- yCols ]+ | i <- [0 .. min ((fst (DX.dimensions df))) n - 1] ]+ sns = TG.snRatioRows (toSN opts) yMatrix+ fes = TG.analyzeSN ad sns+ opts' = TG.optimalLevels fes+ predEta = TG.predictSN fes sns++ printf "Array: %s\n" (T.unpack (OA.oaName (OA.adArray ad)))+ printf "SN type: %s\n" (T.unpack (TG.snTypeName (toSN opts)))+ printf "Inner runs: %d\n" n+ printf "Repetitions per run: %d (columns %s)\n"+ (length yCols) (T.unpack (T.intercalate ", " yCols))+ putStrLn ""++ putStrLn "--- Per-run SN ratios ---"+ mapM_ (\(i, eta) -> printf " Run %2d: SN = %8.3f dB\n" (i :: Int) eta)+ (zip [1..] sns)+ putStrLn ""++ putStrLn "--- Factor effects (mean SN per level) ---"+ mapM_ (printFactorEffect opts') fes+ putStrLn ""++ putStrLn "--- Optimal levels (max SN per factor) ---"+ mapM_ (\(f, lvl, eta) ->+ printf " %-12s = %-12s (SN = %8.3f dB)\n"+ (T.unpack f) (T.unpack (lvText lvl)) eta) opts'+ putStrLn ""+ printf "Predicted SN at optimum (additive model): %.3f dB\n" predEta++ -- ── HTML レポート出力 (--report 指定時) ─────────────────────────────+ case toReport opts of+ Nothing -> return ()+ Just path -> do+ let tr = VTG.TaguchiReport+ { VTG.trTitle = "Taguchi Analysis: "+ <> OA.oaName (OA.adArray ad)+ <> " — "+ <> TG.snTypeName (toSN opts)+ , VTG.trArrayName = OA.oaName (OA.adArray ad)+ , VTG.trSNType = toSN opts+ , VTG.trPerRunSN = sns+ , VTG.trEffects = fes+ , VTG.trOptimal = opts'+ , VTG.trPredicted = predEta+ }+ VTG.renderTaguchiReport path tr+ putStrLn ("Report: " ++ path)+ openInBrowser path+ where+ lvText (OA.LText t) = t+ lvText (OA.LNumeric d)+ | d == fromIntegral (round d :: Integer) = T.pack (show (round d :: Integer))+ | otherwise = T.pack (printf "%g" d)++printFactorEffect :: [(T.Text, OA.LevelValue, Double)] -> TG.FactorEffect -> IO ()+printFactorEffect _opts fe = do+ printf " %s:\n" (T.unpack (TG.feFactor fe))+ let pairs = zip (TG.feLevels fe) (TG.feSNByLevel fe)+ mapM_ (\(lv, eta) ->+ printf " %-12s : %8.3f dB\n"+ (T.unpack (lvShow lv)) eta) pairs+ where+ lvShow (OA.LText t) = t+ lvShow (OA.LNumeric d)+ | d == fromIntegral (round d :: Integer) = T.pack (show (round d :: Integer))+ | otherwise = T.pack (printf "%g" d)++-- ── cross ───────────────────────────────────────────────────────────────++data TgCrossOpts = TgCrossOpts+ { tcInner :: [(T.Text, [T.Text])]+ , tcOuter :: [(T.Text, [T.Text])]+ , tcOut :: Maybe FilePath+ } deriving (Show)++defaultTgCrossOpts :: TgCrossOpts+defaultTgCrossOpts = TgCrossOpts [] [] Nothing++runTaguchiCross :: [String] -> IO ()+runTaguchiCross [] = hPutStrLn stderr "taguchi cross: missing INNER and OUTER array names"+runTaguchiCross [_] = hPutStrLn stderr "taguchi cross: missing OUTER array name"+runTaguchiCross (innerStr : outerStr : rest) =+ case (OA.lookupOA (T.pack innerStr), OA.lookupOA (T.pack outerStr)) of+ (Nothing, _) -> hPutStrLn stderr $+ "taguchi cross: unknown inner array '" ++ innerStr ++ "'"+ (_, Nothing) -> hPutStrLn stderr $+ "taguchi cross: unknown outer array '" ++ outerStr ++ "'"+ (Just innerOA, Just outerOA) ->+ case parseTgCrossOpts rest defaultTgCrossOpts of+ Left err -> hPutStrLn stderr ("taguchi cross: " ++ err)+ Right opts -> doTaguchiCross innerOA outerOA opts++parseTgCrossOpts :: [String] -> TgCrossOpts -> Either String TgCrossOpts+parseTgCrossOpts [] acc = Right acc+parseTgCrossOpts (flag : rest) acc+ | flag `elem` ["-f", "--factor"] = case rest of+ (v : rs) -> case parseFactorSpec v of+ Left err -> Left err+ Right fac -> parseTgCrossOpts rs (acc { tcInner = tcInner acc ++ [fac] })+ [] -> Left "-f/--factor requires NAME=v1,v2,..."+ | flag `elem` ["-fn", "--noise"] = case rest of+ (v : rs) -> case parseFactorSpec v of+ Left err -> Left err+ Right fac -> parseTgCrossOpts rs (acc { tcOuter = tcOuter acc ++ [fac] })+ [] -> Left "-fn/--noise requires NAME=v1,v2,..."+ | flag == "--out" = case rest of+ (v : rs) -> parseTgCrossOpts rs (acc { tcOut = Just v })+ [] -> Left "--out requires a file path"+ | otherwise = Left ("unexpected argument '" ++ flag ++ "'")++doTaguchiCross :: OA.OA -> OA.OA -> TgCrossOpts -> IO ()+doTaguchiCross innerOA outerOA opts = do+ let innerSpecs = [ OA.FactorSpec n (map toLevelValue ls)+ | (n, ls) <- tcInner opts ]+ outerSpecs = [ OA.FactorSpec n (map toLevelValue ls)+ | (n, ls) <- tcOuter opts ]+ case (OA.assignFactors innerOA innerSpecs,+ OA.assignFactors outerOA outerSpecs) of+ (Left err, _) -> hPutStrLn stderr ("inner: " ++ T.unpack err)+ (_, Left err) -> hPutStrLn stderr ("outer: " ++ T.unpack err)+ (Right ai, Right ao) -> do+ let io = TG.makeInnerOuter ai ao+ csv = TG.renderInnerOuterCSV io+ emitText csv (tcOut opts)++-- ---------------------------------------------------------------------------+-- ridge / kernel / spline 共通ヘルパ+-- ---------------------------------------------------------------------------++-- | CSV を読み、x 列(複数可) と y 列(1) を numeric vector で取り出す。+-- 'LoadOpts' を反映 (--no-header / --skip / --comment / --strict)。+loadXY :: LoadOpts -> FilePath -> [T.Text] -> T.Text+ -> IO (Either String (DXD.DataFrame, [V.Vector Double], V.Vector Double))+loadXY lopts path xCols yCol = do+ result <- loadAutoSafeWith lopts path+ case result of+ Left err -> return (Left err)+ Right (df, lg) -> do+ Log.printLogReport lg+ case (mapM (\c -> getDoubleVec c df) xCols, getDoubleVec yCol df) of+ (Just xs, Just y) -> return (Right (df, xs, y))+ _ -> return (Left $ "Numeric column(s) not found: x="+ ++ T.unpack (T.intercalate "," xCols)+ ++ ", y=" ++ T.unpack yCol)++-- | RMSE 計算。+rmseV :: [Double] -> [Double] -> Double+rmseV ys yhat =+ let n = length ys+ sse = sum [ (a - b) ^ (2 :: Int) | (a, b) <- zip ys yhat ]+ in sqrt (sse / fromIntegral (max 1 n))++-- | 散布図 + 滑らか曲線 を出力。+writeSmoothPlot :: OutputFormat -> FilePath -> T.Text+ -> DXD.DataFrame -> T.Text -> T.Text -> SmoothFit -> IO ()+writeSmoothPlot fmt path titleSuffix df xc yc sf =+ scatterWithSmoothFile fmt path+ (defaultConfig (xc <> " vs " <> yc <> " [" <> titleSuffix <> "]"))+ Nothing df xc yc sf++-- | xMin/xMax から評価グリッドを作る。+makeGrid :: V.Vector Double -> Int -> [Double]+makeGrid xs n =+ let lo = V.minimum xs+ hi = V.maximum xs+ in [ lo + fromIntegral i * (hi - lo) / fromIntegral (n - 1)+ | i <- [0 .. n - 1] ]++-- ---------------------------------------------------------------------------+-- ridge subcommand (Ridge / Lasso / Elastic Net)+-- ---------------------------------------------------------------------------++ridgeUsage :: String+ridgeUsage = unlines+ [ "Usage: hanalyze ridge <file> <xcols> <ycol> [options]"+ , ""+ , " <xcols> x column name(s); quote multiple: \"x1 x2\""+ , " <ycol> y column name (single)"+ , ""+ , "Options:"+ , " --penalty TYPE ridge|lasso|elasticnet (default: ridge)"+ , " --lambda L regularization strength (default: 0.1)"+ , " --alpha A ElasticNet L1 mixing in [0,1] (default: 0.5; only with --penalty elasticnet)"+ , " --format FMT html|png|svg (default: html)"+ , " --out FILE scatter+fit output path (default: ridge.html; single x only)"+ , " --report [FILE] build composite HTML report (default: ridge.html)"+ , ""+ , "Examples:"+ , " hanalyze ridge data.csv x y --lambda 0.1"+ , " hanalyze ridge data.csv \"x1 x2 x3\" y --penalty lasso --lambda 0.05"+ , " hanalyze ridge data.csv \"x1 x2\" y --penalty elasticnet --lambda 0.1 --alpha 0.5"+ ]++data RidgeOpts = RidgeOpts+ { roPenalty :: T.Text -- "ridge" / "lasso" / "elasticnet"+ , roLambda :: Double+ , roAlpha :: Double+ , roFormat :: OutputFormat+ , roOut :: FilePath+ , roReport :: Maybe FilePath+ }++defaultRidgeOpts :: RidgeOpts+defaultRidgeOpts = RidgeOpts "ridge" 0.1 0.5 HTML "ridge.html" Nothing++runRidgeCmd :: [String] -> IO ()+runRidgeCmd args0 =+ let (lopts, args) = parseLoadOpts args0+ in case args of+ (file : xColsStr : yColStr : rest) ->+ case parseRidgeOpts rest defaultRidgeOpts of+ Left err -> hPutStrLn stderr ("ridge: " ++ err)+ Right opts -> doRidge file xColsStr yColStr opts lopts+ _ -> putStrLn ridgeUsage++parseRidgeOpts :: [String] -> RidgeOpts -> Either String RidgeOpts+parseRidgeOpts [] acc = Right acc+parseRidgeOpts (flag : rest) acc+ | flag == "--penalty" = case rest of+ (v : rs) | v `elem` ["ridge","lasso","elasticnet"] ->+ parseRidgeOpts rs (acc { roPenalty = T.pack v })+ (v : _) -> Left ("unknown penalty '" ++ v ++ "'")+ [] -> Left "--penalty requires an argument"+ | flag == "--lambda" = case rest of+ (v:rs) -> case reads v :: [(Double, String)] of+ [(d,"")] -> parseRidgeOpts rs (acc { roLambda = d })+ _ -> Left ("invalid --lambda value '" ++ v ++ "'")+ [] -> Left "--lambda requires a value"+ | flag == "--alpha" = case rest of+ (v:rs) -> case reads v :: [(Double, String)] of+ [(d,"")] -> parseRidgeOpts rs (acc { roAlpha = d })+ _ -> Left ("invalid --alpha value '" ++ v ++ "'")+ [] -> Left "--alpha requires a value"+ | flag `elem` ["-f","--format"] = case rest of+ (v:rs) -> case parseFormat v of+ Right f -> parseRidgeOpts rs (acc { roFormat = f })+ Left e -> Left e+ [] -> Left "--format requires an argument"+ | flag == "--out" = case rest of+ (v:rs) -> parseRidgeOpts rs (acc { roOut = v })+ [] -> Left "--out requires a file path"+ | flag == "--report" = case rest of+ (v:rs) | not (null v) && head v /= '-' ->+ parseRidgeOpts rs (acc { roReport = Just v })+ _ -> parseRidgeOpts rest (acc { roReport = Just "ridge.html" })+ | otherwise = Left ("unexpected argument '" ++ flag ++ "'")++doRidge :: FilePath -> String -> String -> RidgeOpts -> LoadOpts -> IO ()+doRidge file xColsStr yColStr opts lopts = do+ let xCols = map T.pack (words xColsStr)+ yCol = T.pack yColStr+ result <- loadXY lopts file xCols yCol+ case result of+ Left err -> hPutStrLn stderr err+ Right (df, xVecs, yVec) -> do+ let n = V.length yVec+ intercept = LA.konst 1 n+ xMat = LA.fromColumns+ (intercept : map (LA.fromList . V.toList) xVecs)+ yLA = LA.fromList (V.toList yVec)+ pen = case roPenalty opts of+ "ridge" -> Reg.L2 (roLambda opts)+ "lasso" -> Reg.L1 (roLambda opts)+ "elasticnet" -> Reg.ElasticNet+ (roLambda opts * roAlpha opts)+ (roLambda opts * (1 - roAlpha opts))+ _ -> Reg.L2 (roLambda opts)+ fit = Reg.fitRegularized pen xMat yLA+ beta = LA.toList (Reg.rfBeta fit)+ yhat = LA.toList (Reg.rfYHat fit)+ ys = V.toList yVec+ rmseVal = rmseV ys yhat+ printf "Loaded %d rows from %s\n" n file+ printf "Penalty: %s, lambda=%g%s\n"+ (T.unpack (roPenalty opts)) (roLambda opts)+ (if roPenalty opts == "elasticnet"+ then ", alpha=" ++ show (roAlpha opts) else "")+ putStrLn ""+ putStrLn "Coefficients:"+ printf " %-30s = %9.4f\n" ("intercept" :: String) (head beta)+ mapM_ (\(i, c, b) ->+ printf " %-30s = %9.4f\n"+ ("β_" ++ show (i :: Int) ++ " (" ++ T.unpack c ++ ")") b)+ (zip3 [1..] xCols (tail beta))+ printf "R² = %.4f\n" (Reg.rfR2 fit)+ printf "|β| > 1e-8: %d / %d (sparsity)\n"+ (Reg.rfNonZero fit) (length beta)+ printf "RMSE (in-sample) = %.4f\n" rmseVal+ -- 単純散布図 + 予測曲線 (1 変数のみ)+ let coeffPairs = zip ("intercept" : xCols)+ (map T.pack (map (printf "%.4f") beta) :: [T.Text])+ coeffNumPairs = zip ("intercept" : xCols) beta+ residuals = LA.toList (Reg.rfResid fit)+ case xCols of+ [xc1] -> do+ let xs = V.toList (head xVecs)+ grid = makeGrid (head xVecs) 100+ gridMat = LA.fromColumns+ [ LA.konst 1 100+ , LA.fromList grid ]+ gridY = LA.toList (Reg.predictRegularized fit gridMat)+ sf = SmoothFit+ { sfX = grid+ , sfFit = gridY+ , sfLower = []+ , sfUpper = []+ , sfHasBand = False+ }+ _ = xs+ _ = coeffPairs+ writeSmoothPlot (roFormat opts) (roOut opts)+ (T.pack ("Regularized: " ++ T.unpack (roPenalty opts)))+ df xc1 yCol sf+ putStrLn ("Plot: " ++ roOut opts)+ openInBrowser (roOut opts)+ -- HTML レポート+ case roReport opts of+ Nothing -> return ()+ Just rpath -> do+ let smooth = RB.SmoothCurve grid gridY [] []+ pathSec = mkRidgePathSection xCols xMat yLA opts+ cfg = ridgeReportConfig opts xCols yCol+ sections =+ [ RB.secDataOverview df xCols yCol+ , RB.secModelOverview (ridgeModelLabel opts)+ (ridgeFormula opts xCols yCol) Nothing+ , RB.secCoefficients coeffNumPairs (Just ("R²", Reg.rfR2 fit))+ , RB.secKeyValue "Fit summary"+ (ridgeFitKVs opts fit beta rmseVal)+ , pathSec+ , RB.secFitScatter xc1 yCol xs ys (Just smooth)+ , RB.secResiduals yhat residuals+ ]+ RB.renderReport rpath cfg sections+ putStrLn ("Report: " ++ rpath)+ openInBrowser rpath+ _ -> do+ putStrLn "(scatter plot skipped for multiple x columns)"+ case roReport opts of+ Nothing -> return ()+ Just rpath -> do+ let pathSec = mkRidgePathSection xCols xMat yLA opts+ cfg = ridgeReportConfig opts xCols yCol+ sections =+ [ RB.secDataOverview df xCols yCol+ , RB.secModelOverview (ridgeModelLabel opts)+ (ridgeFormula opts xCols yCol) Nothing+ , RB.secCoefficients coeffNumPairs (Just ("R²", Reg.rfR2 fit))+ , RB.secKeyValue "Fit summary"+ (ridgeFitKVs opts fit beta rmseVal)+ , pathSec+ , RB.secResiduals yhat residuals+ ]+ RB.renderReport rpath cfg sections+ putStrLn ("Report: " ++ rpath)+ openInBrowser rpath++-- ---------------------------------------------------------------------------+-- kernel subcommand (Nadaraya-Watson / Kernel Ridge / RFF)+-- ---------------------------------------------------------------------------++kernelUsage :: String+kernelUsage = unlines+ [ "Usage: hanalyze kernel <file> <xcol> <ycol> [options]"+ , ""+ , "Options:"+ , " --method M nw|kr|rff (default: kr)"+ , " nw = Nadaraya-Watson"+ , " kr = Kernel Ridge"+ , " rff = Random Fourier Features (RBF)"+ , " --kernel KIND gaussian|epanechnikov|triangular|tricube|uniform"+ , " (default: gaussian; ignored for --method rff)"+ , " --bandwidth H kernel bandwidth h (default: auto via LOO-CV grid)"+ , " --lambda L ridge regularization (default: 0.01; for kr / rff only)"+ , " --features D RFF feature dimension (default: 200; --method rff only)"+ , " --format FMT html|png|svg (default: html)"+ , " --out FILE scatter+fit output path (default: kernel.html)"+ , " --report [FILE] build composite HTML report (default: kernel.html)"+ , ""+ , "Multivariate RFF (--method rff with multiple x columns):"+ , " --group COL group column for color-coded scatter+fit (e.g. name)"+ , " --xaxis COL column to use as horizontal axis in the plot (e.g. t)"+ , " --interactive スライダで副軸を変えると JS が予測曲線を再計算"+ , " (--report と併用、--xaxis の列以外がスライダになる)"+ , " --standardize 入力 X を z-score 化してから fit (スケール差対策)"+ , " --auto-hp HP 自動決定 (default method=loocv)"+ , " --auto-hp-method M loocv (Ridge LOOCV 解析解、推奨) | mlik (周辺尤度最大化)"+ , " (--bandwidth / --lambda は無視される)"+ , ""+ , "Examples:"+ , " hanalyze kernel data.csv x y --method kr --bandwidth 0.5"+ , " hanalyze kernel data.csv x y --method nw # auto-bandwidth via LOO-CV"+ , " hanalyze kernel data.csv x y --method rff --features 200"+ , " # 多変量 RFF (melted データに対して):"+ , " hanalyze kernel data/io/melted_sample.csv \"x1 t\" y --method rff \\"+ , " --features 200 --bandwidth 1.0 --lambda 0.001 \\"+ , " --group name --xaxis t --out plot.html"+ ]++data KernelOpts = KernelOpts+ { koMethod :: T.Text -- "nw" / "kr" / "rff"+ , koKernel :: Kern.Kernel -- Gaussian / Epanechnikov / ...+ , koBandwidth :: Maybe Double+ , koLambda :: Double+ , koFeatures :: Int+ , koFormat :: OutputFormat+ , koOut :: FilePath+ , koReport :: Maybe FilePath+ , koGroup :: Maybe T.Text -- 多変量 RFF プロット用 group 列+ , koXAxis :: Maybe T.Text -- 多変量 RFF プロット用 横軸列名+ , koInteractive :: Bool -- インタラクティブ予測 (--report と併用)+ , koStandardize :: Bool -- 入力標準化 (Phase 4)+ , koAutoHP :: Bool -- HP 自動決定+ , koAutoHPMethod :: T.Text -- "loocv" / "mlik" (default loocv = 速い)+ }++defaultKernelOpts :: KernelOpts+defaultKernelOpts = KernelOpts+ { koMethod = "kr"+ , koKernel = Kern.Gaussian+ , koBandwidth = Nothing+ , koLambda = 0.01+ , koFeatures = 200+ , koFormat = HTML+ , koOut = "kernel.html"+ , koReport = Nothing+ , koGroup = Nothing+ , koXAxis = Nothing+ , koInteractive = False+ , koStandardize = False+ , koAutoHP = False+ , koAutoHPMethod = "loocv"+ }++parseKernelKind :: String -> Either String Kern.Kernel+parseKernelKind s = case s of+ "gaussian" -> Right Kern.Gaussian+ "epanechnikov" -> Right Kern.Epanechnikov+ "triangular" -> Right Kern.Triangular+ "tricube" -> Right Kern.TriCube+ "uniform" -> Right Kern.Uniform+ _ -> Left ("unknown kernel '" ++ s ++ "'")++runKernelCmd :: [String] -> IO ()+runKernelCmd args0 =+ let (lopts, args) = parseLoadOpts args0+ in case args of+ (file : xColStr : yColStr : rest) ->+ case parseKernelOpts rest defaultKernelOpts of+ Left err -> hPutStrLn stderr ("kernel: " ++ err)+ Right opts -> doKernel file xColStr yColStr opts lopts+ _ -> putStrLn kernelUsage++parseKernelOpts :: [String] -> KernelOpts -> Either String KernelOpts+parseKernelOpts [] acc = Right acc+parseKernelOpts (flag : rest) acc+ | flag == "--method" = case rest of+ (v:rs) | v `elem` ["nw","kr","rff"] ->+ parseKernelOpts rs (acc { koMethod = T.pack v })+ (v:_) -> Left ("unknown method '" ++ v ++ "'")+ [] -> Left "--method requires an argument"+ | flag == "--kernel" = case rest of+ (v:rs) -> case parseKernelKind v of+ Right k -> parseKernelOpts rs (acc { koKernel = k })+ Left e -> Left e+ [] -> Left "--kernel requires an argument"+ | flag == "--bandwidth" = case rest of+ (v:rs) -> case reads v :: [(Double, String)] of+ [(d,"")] -> parseKernelOpts rs (acc { koBandwidth = Just d })+ _ -> Left ("invalid --bandwidth '" ++ v ++ "'")+ [] -> Left "--bandwidth requires a value"+ | flag == "--lambda" = case rest of+ (v:rs) -> case reads v :: [(Double, String)] of+ [(d,"")] -> parseKernelOpts rs (acc { koLambda = d })+ _ -> Left ("invalid --lambda '" ++ v ++ "'")+ [] -> Left "--lambda requires a value"+ | flag == "--features" = case rest of+ (v:rs) -> case reads v :: [(Int, String)] of+ [(d,"")] -> parseKernelOpts rs (acc { koFeatures = d })+ _ -> Left ("invalid --features '" ++ v ++ "'")+ [] -> Left "--features requires a value"+ | flag `elem` ["-f","--format"] = case rest of+ (v:rs) -> case parseFormat v of+ Right f -> parseKernelOpts rs (acc { koFormat = f })+ Left e -> Left e+ [] -> Left "--format requires an argument"+ | flag == "--out" = case rest of+ (v:rs) -> parseKernelOpts rs (acc { koOut = v })+ [] -> Left "--out requires a file path"+ | flag == "--report" = case rest of+ (v:rs) | not (null v) && head v /= '-' ->+ parseKernelOpts rs (acc { koReport = Just v })+ _ -> parseKernelOpts rest (acc { koReport = Just "kernel.html" })+ | flag == "--group" = case rest of+ (v:rs) -> parseKernelOpts rs (acc { koGroup = Just (T.pack v) })+ [] -> Left "--group requires a column name"+ | flag == "--xaxis" = case rest of+ (v:rs) -> parseKernelOpts rs (acc { koXAxis = Just (T.pack v) })+ [] -> Left "--xaxis requires a column name"+ | flag == "--interactive" =+ parseKernelOpts rest (acc { koInteractive = True })+ | flag == "--standardize" =+ parseKernelOpts rest (acc { koStandardize = True })+ | flag == "--auto-hp" =+ parseKernelOpts rest (acc { koAutoHP = True })+ | flag == "--auto-hp-method" = case rest of+ (v:rs) | v `elem` ["loocv", "mlik"] ->+ parseKernelOpts rs (acc { koAutoHP = True+ , koAutoHPMethod = T.pack v })+ (v:_) -> Left ("unknown --auto-hp-method '" ++ v ++ "' (choose loocv|mlik)")+ [] -> Left "--auto-hp-method requires loocv|mlik"+ | otherwise = Left ("unexpected argument '" ++ flag ++ "'")++doKernel :: FilePath -> String -> String -> KernelOpts -> LoadOpts -> IO ()+doKernel file xColStr yColStr opts lopts = do+ let xCols = map T.pack (words xColStr)+ yCol = T.pack yColStr+ case xCols of+ [] -> hPutStrLn stderr "kernel: x 列が指定されていません"+ [xCol] -> do+ result <- loadXY lopts file [xCol] yCol+ case result of+ Left err -> hPutStrLn stderr err+ Right (df, [xVec], yVec) ->+ runKernelOn df xCol yCol xVec yVec opts+ Right _ -> hPutStrLn stderr "kernel: expected single x column"+ _multiple -> case koMethod opts of+ "rff" -> do+ result <- loadXY lopts file xCols yCol+ case result of+ Left err -> hPutStrLn stderr err+ Right (df, xVecs, yVec) ->+ runKernelMV df xCols yCol xVecs yVec opts+ "kr" -> do+ result <- loadXY lopts file xCols yCol+ case result of+ Left err -> hPutStrLn stderr err+ Right (_, xVecs, yVec) ->+ runKernelMVKR xCols yCol xVecs yVec opts+ "nw" -> do+ result <- loadXY lopts file xCols yCol+ case result of+ Left err -> hPutStrLn stderr err+ Right (_, xVecs, yVec) ->+ runKernelMVNW xCols yCol xVecs yVec opts+ m -> hPutStrLn stderr $+ "kernel --method " ++ T.unpack m+ ++ " (unknown method)"++-- | Multi-input Kernel Ridge (Phase K5) — fit and report training metrics.+-- 多次元 X (n×p) を取り、Hanalyze.Model.KernelRegression.kernelRidgeMV で fit。+-- 予測図は生成しない (多次元のため)、R² と RMSE をログ出力。+runKernelMVKR+ :: [T.Text] -> T.Text -> [V.Vector Double] -> V.Vector Double+ -> KernelOpts -> IO ()+runKernelMVKR xCols _yCol xVecs yVec opts = do+ let n = V.length yVec+ p = length xCols+ xMat = LA.fromColumns (map (LA.fromList . V.toList) xVecs)+ yMat = LA.asColumn (LA.fromList (V.toList yVec))+ ker = koKernel opts+ h = case koBandwidth opts of+ Just b -> b+ Nothing -> 1.0+ lam = koLambda opts+ fit = Kern.kernelRidgeMV ker h lam xMat yMat+ yhat = Kern.fittedKernelRidgeMV fit+ ss = LA.sumElements ((yMat - yhat) ** 2)+ muY = LA.sumElements yMat / fromIntegral n+ stTot = LA.sumElements ((yMat - LA.konst muY (n, 1)) ** 2)+ r2 = 1 - ss / stTot+ rmse = sqrt (ss / fromIntegral n)+ printf "Loaded %d rows × %d features (%s); method=kr (multivariate)\n"+ n p (T.unpack (T.intercalate "," xCols))+ printf " bandwidth h = %.4g, lambda = %.4g, kernel = %s\n"+ h lam (show ker)+ printf " R² (train) = %.4f\n" r2+ printf " RMSE (train) = %.4f\n" rmse++-- | Multi-input Nadaraya-Watson (Phase K5) — fit and report training metrics.+runKernelMVNW+ :: [T.Text] -> T.Text -> [V.Vector Double] -> V.Vector Double+ -> KernelOpts -> IO ()+runKernelMVNW xCols _yCol xVecs yVec opts = do+ let n = V.length yVec+ p = length xCols+ xMat = LA.fromColumns (map (LA.fromList . V.toList) xVecs)+ yMat = LA.asColumn (LA.fromList (V.toList yVec))+ ker = koKernel opts+ h = case koBandwidth opts of+ Just b -> b+ Nothing -> 1.0+ yhat = Kern.nwRegressionMV ker h xMat yMat xMat+ ss = LA.sumElements ((yMat - yhat) ** 2)+ muY = LA.sumElements yMat / fromIntegral n+ stTot = LA.sumElements ((yMat - LA.konst muY (n, 1)) ** 2)+ r2 = 1 - ss / stTot+ rmse = sqrt (ss / fromIntegral n)+ printf "Loaded %d rows × %d features (%s); method=nw (multivariate)\n"+ n p (T.unpack (T.intercalate "," xCols))+ printf " bandwidth h = %.4g, kernel = %s\n" h (show ker)+ printf " R² (train) = %.4f\n" r2+ printf " RMSE (train) = %.4f\n" rmse++-- | 多変量 RFF Ridge を走らせる (Phase B-RFF)。+-- '--group' / '--xaxis' が指定されていれば、グループ別観測点 + 予測曲線の+-- 散布図を出力する。+-- '--standardize' / '--auto-hp' で前処理 / HP 自動決定。+runKernelMV+ :: DXD.DataFrame -> [T.Text] -> T.Text+ -> [V.Vector Double] -> V.Vector Double+ -> KernelOpts -> IO ()+runKernelMV df xCols yCol xVecs yVec opts = do+ let n = V.length yVec+ p = length xCols+ cols = map V.toList xVecs+ xMatRaw = LA.fromColumns (map LA.fromList cols)+ ys = V.toList yVec+ yV = LA.fromList ys+ printf "Loaded %d rows × %d features (%s); method=rff (multivariate)\n"+ n p (T.unpack (T.intercalate "," xCols))++ -- ステップ 1: 標準化 (タイマー付き)+ (tStd, (stdr, xMat)) <- timed $ do+ let s = if koStandardize opts+ then Std.fitStandardizer xMatRaw+ else Std.identityStandardizer p+ xm = if koStandardize opts+ then Std.applyStandardizer s xMatRaw+ else xMatRaw+ return (s, xm)+ if koStandardize opts+ then do+ putStrLn " Standardize: ON"+ printf " μ = [%s]\n" (T.unpack (T.intercalate ", " (map NF.fmtNumT (Std.stMu stdr))))+ printf " σ = [%s]\n" (T.unpack (T.intercalate ", " (map NF.fmtNumT (Std.stSd stdr))))+ else putStrLn " Standardize: OFF"++ -- ステップ 2: HP の決定 (タイマー付き)+ hpGen <- createSystemRandom+ (tHP, (ell, lam, sigF)) <- timed $+ if koAutoHP opts+ then case koAutoHPMethod opts of+ "loocv" -> do+ putStrLn " Auto-HP (LOOCV): RFF Ridge の解析的 LOO を最小化中..."+ res <- RFF.gridSearchLOOCVRBFMV p (koFeatures opts) xMat yV Nothing hpGen+ let ellOpt = RFF.lcEll res+ sfOpt = RFF.lcSigmaF res+ lamOpt = RFF.lcLambda res+ ellOpt `seq` sfOpt `seq` lamOpt `seq` return ()+ printf " ℓ = %s\n" (NF.fmtNum ellOpt)+ printf " σ_f = %s\n" (NF.fmtNum sfOpt)+ printf " λ = %s\n" (NF.fmtNum lamOpt)+ printf " LOOCV = %s (グリッド %d 点評価)\n"+ (NF.fmtNum (RFF.lcLOOCV res)) (RFF.lcGridPts res)+ return (ellOpt, lamOpt, sfOpt)+ _ -> do -- "mlik"+ putStrLn " Auto-HP (周辺尤度): Cholesky で marg-lik を最大化中..."+ let res = RFF.maximizeMarginalLikRBFMV xMat yV Nothing+ ellOpt = RFF.mlEll res+ sfOpt = RFF.mlSigmaF res+ snOpt = RFF.mlSigmaN res+ lamOpt = snOpt * snOpt+ ellOpt `seq` sfOpt `seq` snOpt `seq` return ()+ printf " ℓ = %s\n" (NF.fmtNum ellOpt)+ printf " σ_f = %s\n" (NF.fmtNum sfOpt)+ printf " σ_n = %s (λ = σ_n² = %s)\n"+ (NF.fmtNum snOpt) (NF.fmtNum lamOpt)+ printf " log_mlik = %s (グリッド %d 点評価)\n"+ (NF.fmtNum (RFF.mlLogMlik res)) (RFF.mlGridPts res)+ return (ellOpt, lamOpt, sfOpt)+ else do+ let ell0 = case koBandwidth opts of+ Just h -> h+ Nothing -> defaultLengthScale (map LA.toList (LA.toColumns xMat))+ printf " ell=%s lambda=%s\n"+ (NF.fmtNum ell0) (NF.fmtNum (koLambda opts))+ return (ell0, koLambda opts, 1.0)++ let d = koFeatures opts+ printf " D=%d\n" d++ -- ステップ 3: RFF サンプリング + Ridge fit + 評価 (各タイマー付き)+ gen <- createSystemRandom+ (tSample, feats) <- timed (RFF.sampleRFFRBFMV p d ell sigF gen)+ (tFit, fit) <- timed $ do+ let f = RFF.rffRidgeMV feats xMat ys lam+ LA.size (RFF.rffrmvWeights f) `seq` return f+ let yhat = RFF.predictRFFRidgeMV fit xMat+ sse = sum (zipWith (\a b -> (a - b)^(2::Int)) ys yhat)+ sst = let m = sum ys / fromIntegral (max 1 (length ys))+ in sum [(y - m)^(2::Int) | y <- ys]+ r2 = if sst < 1e-12 then 0 else 1 - sse / sst+ printf "RFF (multivariate) Ridge fit:\n"+ printf " R^2 = %s\n" (NF.fmtNum r2)+ printf " RMSE = %s\n" (NF.fmtNum (sqrt (sse / fromIntegral n)))++ putStrLn ""+ putStrLn "Profiling (cumulative wall time):"+ printPhase "Standardize" tStd+ printPhase "Auto-HP" tHP+ printPhase "Sample RFF" tSample+ printPhase "Fit Ridge" tFit++ -- --group + --xaxis が両方指定されていればプロット+ case (koGroup opts, koXAxis opts) of+ (Just gCol, Just xCol) -> do+ let outPath = koOut opts+ fmt = koFormat opts+ (tPlot, _) <- timed (writeMVPlot fmt outPath df gCol xCol xCols yCol fit stdr cols ys)+ putStrLn $ "Plot: " ++ outPath+ printPhase "Plot" tPlot+ -- --report 指定時は ReportBuilder で統合 HTML を出力+ case koReport opts of+ Just rpath -> do+ (tRep, _) <- timed $ do+ let rep = RI.RFFMVReport+ { RI.rfmvFit = fit+ , RI.rfmvGroup = gCol+ , RI.rfmvXAxis = xCol+ , RI.rfmvInteractive = koInteractive opts+ , RI.rfmvStandardizer =+ if koStandardize opts then Just stdr else Nothing+ }+ cfg = RB.defaultReportConfig+ (yCol <> " — Multivariate RFF Ridge"+ <> if koInteractive opts then " (interactive)" else "")+ secs = RB.toReport cfg df xCols yCol rep+ RB.renderReport rpath cfg secs+ putStrLn $ "Report: " ++ rpath+ printPhase "Render report" tRep+ Nothing -> return ()+ _ -> putStrLn+ "Plot skipped (use --group COL --xaxis COL to draw scatter+fit by group)"++-- | name (group) ごとに観測点と予測曲線をプロット。+-- 標準化 ON のときは予測グリッドを raw → 標準化空間に変換してから predict。+-- 横軸 / 観測点は raw 単位で表示する。+writeMVPlot+ :: OutputFormat -> FilePath+ -> DXD.DataFrame+ -> T.Text -> T.Text -> [T.Text] -> T.Text+ -> RFF.RFFRidgeFitMV+ -> Std.Standardizer+ -> [[Double]] -- ^ raw cols+ -> [Double]+ -> IO ()+writeMVPlot fmt path df gCol xCol xCols yCol fit stdr cols ys = do+ case getMaybeTextVec gCol df of+ Nothing -> hPutStrLn stderr $+ "plot: group column '" ++ T.unpack gCol ++ "' not found"+ Just gv ->+ let groups = [ maybe "" id g | g <- V.toList gv ]+ xColIdx = case [ i | (i, c) <- zip [0..] xCols, c == xCol ] of+ (i:_) -> i+ [] -> 0+ xValuesAll = cols !! xColIdx+ xMin = minimum xValuesAll+ xMax = maximum xValuesAll+ ngrid = 100+ xGrid = [ xMin + fromIntegral i * (xMax - xMin) / fromIntegral (ngrid - 1)+ | i <- [0 .. ngrid - 1] ]+ ptData = zip3 groups xValuesAll ys+ uniqGroups = uniq groups+ rowsForGroup g = [ i | (i, gg) <- zip [0..] groups, gg == g ]+ repValues g = [ (cols !! j) !! head (rowsForGroup g)+ | j <- [0 .. length xCols - 1] ]+ mkLineData g =+ let rep = repValues g+ -- raw 値で row を組む+ makeRowRaw t =+ [ if j == xColIdx then t else rep !! j+ | j <- [0 .. length xCols - 1] ]+ xMatRaw = LA.fromLists [ makeRowRaw t | t <- xGrid ]+ -- 標準化空間に変換してから predict+ xMatStd = Std.applyStandardizer stdr xMatRaw+ ys' = RFF.predictRFFRidgeMV fit xMatStd+ in [ (g, t, y') | (t, y') <- zip xGrid ys' ]+ lnData = concatMap mkLineData uniqGroups+ plotCfg = (defaultConfig (yCol <> " by " <> gCol))+ { plotWidth = 720, plotHeight = 480 }+ in scatterWithGroupsFile fmt path plotCfg xCol yCol ptData lnData++uniq :: Ord a => [a] -> [a]+uniq [] = []+uniq (x:xs) = x : uniq (filter (/= x) xs)++-- | 各列の標準偏差の幾何平均で長さスケールを推定 (median heuristic 簡易版)。+defaultLengthScale :: [[Double]] -> Double+defaultLengthScale cols =+ let stds = [ std c | c <- cols, length c > 1 ]+ std xs = let n = fromIntegral (length xs)+ m = sum xs / n+ v = sum [ (x - m)^(2::Int) | x <- xs ] / max 1 (n - 1)+ in sqrt v+ g = product stds ** (1.0 / fromIntegral (max 1 (length stds)))+ in if g <= 0 then 1.0 else g++runKernelOn :: DXD.DataFrame -> T.Text -> T.Text -> V.Vector Double -> V.Vector Double+ -> KernelOpts -> IO ()+runKernelOn df xCol yCol xVec yVec opts = do+ let n = V.length xVec+ method = koMethod opts+ ker = koKernel opts+ grid = makeGrid xVec 100+ gridV = V.fromList grid+ printf "Loaded %d rows; method=%s, kernel=%s\n"+ n (T.unpack method) (show ker)++ -- Bandwidth selection+ h <- case koBandwidth opts of+ Just hVal -> do+ printf "Bandwidth (specified): h = %.4f\n" hVal+ return hVal+ Nothing -> do+ let xMin = V.minimum xVec+ xMax = V.maximum xVec+ range = xMax - xMin+ hCands = [range/40, range/20, range/10, range/5, range/2.5]+ (bestH, bestRMSE) = Kern.gridSearchBandwidth ker xVec yVec hCands+ printf "Bandwidth (LOO-CV best): h = %.4f (CV-RMSE = %.4f)\n"+ bestH bestRMSE+ return bestH++ -- Fit + predict on grid+ (gridY, sumStr) <- case method of+ "nw" -> do+ let ys = Kern.nwRegression ker h xVec yVec gridV+ return (V.toList ys, "Nadaraya-Watson, h=" ++ show h)+ "kr" -> do+ let lam = koLambda opts+ fit = Kern.kernelRidge ker h lam xVec yVec+ ys = Kern.predictKernelRidge fit gridV+ return (V.toList ys+ , "Kernel Ridge, h=" ++ show h ++ ", lambda=" ++ show lam)+ "rff" -> do+ gen <- createSystemRandom+ feats <- RFF.sampleRFFRBF (koFeatures opts) h 1.0 gen+ let lam = koLambda opts+ fit = RFF.rffRidge feats (V.toList xVec) (V.toList yVec) lam+ ys = RFF.predictRFFRidge fit grid+ return (ys, "RFF, D=" ++ show (koFeatures opts)+ ++ ", h=" ++ show h ++ ", lambda=" ++ show lam)+ _ -> error "unreachable"++ -- In-sample RMSE+ let predictX :: V.Vector Double -> [Double]+ predictX xs = case method of+ "nw" -> V.toList (Kern.nwRegression ker h xVec yVec xs)+ "kr" -> V.toList (Kern.predictKernelRidge+ (Kern.kernelRidge ker h (koLambda opts) xVec yVec) xs)+ _ -> [] -- rff requires gen; skip in-sample for now+ ys = V.toList yVec+ case method of+ "rff" -> printf "Predictions on %d test points; in-sample RMSE skipped (RFF re-samples)\n"+ (length grid)+ _ -> printf "RMSE (in-sample) = %.4f\n" (rmseV ys (predictX xVec))+ putStrLn $ "(" ++ sumStr ++ ")"++ -- Plot+ let sf = SmoothFit+ { sfX = grid+ , sfFit = gridY+ , sfLower = []+ , sfUpper = []+ , sfHasBand = False+ }+ writeSmoothPlot (koFormat opts) (koOut opts)+ (T.pack ("Kernel: " ++ T.unpack method)) df xCol yCol sf+ putStrLn ("Plot: " ++ koOut opts)+ openInBrowser (koOut opts)++ -- HTML レポート (--report)+ case koReport opts of+ Nothing -> return ()+ Just rpath -> do+ let xs = V.toList xVec+ ys = V.toList yVec+ smooth = RB.SmoothCurve grid gridY [] []+ modelLbl = "Kernel regression (" <> method <> ")"+ formula = T.pack (T.unpack yCol ++ " ~ f(" ++ T.unpack xCol ++ ")")+ cfg = RB.defaultReportConfig+ ("Kernel regression — " <> yCol <> " ~ " <> xCol)+ baseKVs =+ [ ("Method", method)+ , ("Kernel", T.pack (show ker))+ , ("Bandwidth", T.pack (printf "%.4f" h))+ ]+ extraKVs = case method of+ "kr" -> [("Lambda", T.pack (printf "%g" (koLambda opts)))]+ "rff" -> [("Features", T.pack (show (koFeatures opts)))+ ,("Lambda", T.pack (printf "%g" (koLambda opts)))]+ _ -> []+ sections =+ [ RB.secDataOverview df [xCol] yCol+ , RB.secModelOverview modelLbl formula Nothing+ , RB.secKeyValue "Fit summary" (baseKVs ++ extraKVs)+ , RB.secFitScatter xCol yCol xs ys (Just smooth)+ ]+ RB.renderReport rpath cfg sections+ putStrLn ("Report: " ++ rpath)+ openInBrowser rpath++-- ---------------------------------------------------------------------------+-- spline subcommand+-- ---------------------------------------------------------------------------++splineUsage :: String+splineUsage = unlines+ [ "Usage: hanalyze spline <file> <xcol> <ycol> [options]"+ , ""+ , "Options:"+ , " --type T bspline|natural (default: bspline)"+ , " --knots N number of internal knots (default: 5)"+ , " --degree D B-spline degree (default: 3 = cubic)"+ , " --format FMT html|png|svg (default: html)"+ , " --out FILE scatter+fit output path (default: spline.html)"+ , " --report [FILE] build composite HTML report (default: spline.html)"+ , ""+ , "Examples:"+ , " hanalyze spline data.csv x y --knots 8"+ , " hanalyze spline data.csv x y --type natural"+ , " hanalyze spline data.csv x y --type bspline --degree 3 --knots 10"+ ]++data SplineOpts = SplineOpts+ { soType :: T.Text+ , soKnots :: Int+ , soDegree :: Int+ , soFormat :: OutputFormat+ , soOut :: FilePath+ , soReport :: Maybe FilePath+ }++defaultSplineOpts :: SplineOpts+defaultSplineOpts = SplineOpts "bspline" 5 3 HTML "spline.html" Nothing++runSplineCmd :: [String] -> IO ()+runSplineCmd args0 =+ let (lopts, args) = parseLoadOpts args0+ in case args of+ (file : xColStr : yColStr : rest) ->+ case parseSplineOpts rest defaultSplineOpts of+ Left err -> hPutStrLn stderr ("spline: " ++ err)+ Right opts -> doSpline file xColStr yColStr opts lopts+ _ -> putStrLn splineUsage++parseSplineOpts :: [String] -> SplineOpts -> Either String SplineOpts+parseSplineOpts [] acc = Right acc+parseSplineOpts (flag : rest) acc+ | flag == "--type" = case rest of+ (v:rs) | v `elem` ["bspline","natural"] ->+ parseSplineOpts rs (acc { soType = T.pack v })+ (v:_) -> Left ("unknown spline type '" ++ v ++ "'")+ [] -> Left "--type requires an argument"+ | flag == "--knots" = case rest of+ (v:rs) -> case reads v :: [(Int, String)] of+ [(d,"")] -> parseSplineOpts rs (acc { soKnots = d })+ _ -> Left ("invalid --knots '" ++ v ++ "'")+ [] -> Left "--knots requires a value"+ | flag == "--degree" = case rest of+ (v:rs) -> case reads v :: [(Int, String)] of+ [(d,"")] -> parseSplineOpts rs (acc { soDegree = d })+ _ -> Left ("invalid --degree '" ++ v ++ "'")+ [] -> Left "--degree requires a value"+ | flag `elem` ["-f","--format"] = case rest of+ (v:rs) -> case parseFormat v of+ Right f -> parseSplineOpts rs (acc { soFormat = f })+ Left e -> Left e+ [] -> Left "--format requires an argument"+ | flag == "--out" = case rest of+ (v:rs) -> parseSplineOpts rs (acc { soOut = v })+ [] -> Left "--out requires a file path"+ | flag == "--report" = case rest of+ (v:rs) | not (null v) && head v /= '-' ->+ parseSplineOpts rs (acc { soReport = Just v })+ _ -> parseSplineOpts rest (acc { soReport = Just "spline.html" })+ | otherwise = Left ("unexpected argument '" ++ flag ++ "'")++doSpline :: FilePath -> String -> String -> SplineOpts -> LoadOpts -> IO ()+doSpline file xColStr yColStr opts lopts = do+ let xCol = T.pack xColStr+ yCol = T.pack yColStr+ result <- loadXY lopts file [xCol] yCol+ case result of+ Left err -> hPutStrLn stderr err+ Right (df, [xVec], yVec) -> do+ let kind = case soType opts of+ "natural" -> Spl.NaturalCubic+ _ -> Spl.BSpline (soDegree opts)+ k = soKnots opts+ xMin = V.minimum xVec+ xMax = V.maximum xVec+ knots = [ xMin + fromIntegral i * (xMax - xMin) / fromIntegral (k + 1)+ | i <- [1 .. k] ]+ fit = Spl.fitSpline kind knots xVec yVec+ grid = makeGrid xVec 100+ gridV = V.fromList grid+ gridY = V.toList (Spl.predictSpline fit gridV)+ n = V.length xVec+ ys = V.toList yVec+ yhatIn = V.toList (Spl.predictSpline fit xVec)+ rmseVal = rmseV ys yhatIn+ printf "Loaded %d rows; type=%s, knots=%d%s\n"+ n (T.unpack (soType opts)) k+ (if soType opts == "bspline"+ then ", degree=" ++ show (soDegree opts) else "")+ printf "RMSE (in-sample) = %.4f\n" rmseVal+ let sf = SmoothFit+ { sfX = grid+ , sfFit = gridY+ , sfLower = []+ , sfUpper = []+ , sfHasBand = False+ }+ writeSmoothPlot (soFormat opts) (soOut opts)+ (T.pack ("Spline: " ++ T.unpack (soType opts))) df xCol yCol sf+ putStrLn ("Plot: " ++ soOut opts)+ openInBrowser (soOut opts)+ -- HTML レポート (--report)+ case soReport opts of+ Nothing -> return ()+ Just rpath -> do+ let smooth = RB.SmoothCurve grid gridY [] []+ modelLbl = "Spline regression (" <> soType opts <> ")"+ formula = T.pack (T.unpack yCol ++ " ~ s("+ ++ T.unpack xCol ++ "; knots="+ ++ show k ++ ")")+ cfg = RB.defaultReportConfig+ ("Spline regression — " <> yCol <> " ~ " <> xCol)+ sections =+ [ RB.secDataOverview df [xCol] yCol+ , RB.secModelOverview modelLbl formula Nothing+ , RB.secKeyValue "Fit summary"+ [ ("Type", soType opts)+ , ("Knots", T.pack (show k))+ , ("Degree", T.pack (show (soDegree opts)))+ , ("RMSE (in-sample)", T.pack (printf "%.4f" rmseVal))+ ]+ , RB.secFitScatter xCol yCol (V.toList xVec) ys+ (Just smooth)+ , RB.secResiduals yhatIn (zipWith (-) ys yhatIn)+ ]+ RB.renderReport rpath cfg sections+ putStrLn ("Report: " ++ rpath)+ openInBrowser rpath+ Right _ -> hPutStrLn stderr "spline: expected single x column"++-- ---------------------------------------------------------------------------+-- ridge report ヘルパ+-- ---------------------------------------------------------------------------++ridgeModelLabel :: RidgeOpts -> T.Text+ridgeModelLabel opts =+ "Regularized regression (" <> roPenalty opts <> ")"++ridgeFormula :: RidgeOpts -> [T.Text] -> T.Text -> T.Text+ridgeFormula opts xCols yCol =+ T.pack (T.unpack yCol ++ " ~ "+ ++ intercalate " + " (map T.unpack xCols)+ ++ " (lambda=" ++ show (roLambda opts) ++ ")")++ridgeReportConfig :: RidgeOpts -> [T.Text] -> T.Text -> RB.ReportConfig+ridgeReportConfig _opts xCols yCol = RB.defaultReportConfig+ ("Regularized regression — "+ <> yCol <> " ~ " <> T.intercalate " + " xCols)++ridgeFitKVs :: RidgeOpts -> Reg.RegFit -> [Double] -> Double -> [(T.Text, T.Text)]+ridgeFitKVs opts fit beta rmseVal =+ [ ("RMSE (in-sample)", T.pack (printf "%.4f" rmseVal))+ , ("|β| > 1e-8", T.pack (show (Reg.rfNonZero fit) <> " / "+ <> show (length beta)))+ , ("Penalty", roPenalty opts)+ , ("Lambda", T.pack (printf "%g" (roLambda opts)))+ ]++-- | Regularization path: λ を 1e-4 .. 1e2 で対数スケール掃引、+-- 各 λ で fit して係数を集める。intercept は除外して可視化。+mkRidgePathSection :: [T.Text] -> LA.Matrix Double -> LA.Vector Double+ -> RidgeOpts -> RB.ReportSection+mkRidgePathSection xCols xMat yLA opts =+ let lambdas = [10 ** (-4 + 0.1 * fromIntegral i) | i <- [0 .. 60 :: Int]]+ mkPen lam = case roPenalty opts of+ "ridge" -> Reg.L2 lam+ "lasso" -> Reg.L1 lam+ "elasticnet" -> Reg.ElasticNet (lam * roAlpha opts)+ (lam * (1 - roAlpha opts))+ _ -> Reg.L2 lam+ path = Reg.regularizationPath mkPen lambdas xMat yLA+ -- intercept (係数 0) を除外+ pathNoInt = [ (lam, drop 1 coefs) | (lam, coefs) <- path ]+ title = "Regularization path (" <> roPenalty opts <> ")"+ spec = RB.regPathSpec xCols pathNoInt+ in RB.secVega title spec++-- ---------------------------------------------------------------------------+-- quantile subcommand+-- ---------------------------------------------------------------------------++quantileUsage :: String+quantileUsage = unlines+ [ "Usage: hanalyze quantile <file> <xcols> <ycol> [options]"+ , ""+ , " <xcols> x column name(s); quote multiple: \"x1 x2\""+ , " <ycol> y column name (single)"+ , ""+ , "Options:"+ , " --tau T quantile in (0, 1) (default: 0.5 = median)"+ , " --taus T1,T2,... overlay multiple quantiles in the report (e.g. 0.1,0.5,0.9)"+ , " --format FMT html|png|svg (default: html)"+ , " --out FILE scatter+fit output path (default: quantile.html)"+ , " --report [FILE] build composite HTML report (default: quantile.html)"+ , ""+ , "Examples:"+ , " hanalyze quantile data.csv x y --tau 0.5"+ , " hanalyze quantile data.csv x y --taus 0.1,0.5,0.9 --report"+ ]++data QuantileOpts = QuantileOpts+ { qoTau :: Double+ , qoTaus :: [Double] -- when not empty, overlay multiple quantiles+ , qoFormat :: OutputFormat+ , qoOut :: FilePath+ , qoReport :: Maybe FilePath+ }++defaultQuantileOpts :: QuantileOpts+defaultQuantileOpts = QuantileOpts 0.5 [] HTML "quantile.html" Nothing++runQuantileCmd :: [String] -> IO ()+runQuantileCmd args0 =+ let (lopts, args) = parseLoadOpts args0+ in case args of+ (file : xColsStr : yColStr : rest) ->+ case parseQuantileOpts rest defaultQuantileOpts of+ Left err -> hPutStrLn stderr ("quantile: " ++ err)+ Right opts -> doQuantile file xColsStr yColStr opts lopts+ _ -> putStrLn quantileUsage++parseQuantileOpts :: [String] -> QuantileOpts -> Either String QuantileOpts+parseQuantileOpts [] acc = Right acc+parseQuantileOpts (flag : rest) acc+ | flag == "--tau" = case rest of+ (v:rs) -> case reads v :: [(Double, String)] of+ [(d,"")] | d > 0, d < 1 -> parseQuantileOpts rs (acc { qoTau = d })+ _ -> Left ("invalid --tau '" ++ v ++ "' (must be in (0,1))")+ [] -> Left "--tau requires a value"+ | flag == "--taus" = case rest of+ (v:rs) ->+ let parts = filter (not . null) (splitOnComma v)+ in case mapM (\s -> case reads s :: [(Double, String)] of+ [(d,"")] | d > 0, d < 1 -> Just d+ _ -> Nothing) parts of+ Just ds -> parseQuantileOpts rs (acc { qoTaus = ds })+ Nothing -> Left ("invalid --taus '" ++ v+ ++ "' (comma-separated values in (0,1))")+ [] -> Left "--taus requires a value"+ | flag `elem` ["-f","--format"] = case rest of+ (v:rs) -> case parseFormat v of+ Right f -> parseQuantileOpts rs (acc { qoFormat = f })+ Left e -> Left e+ [] -> Left "--format requires an argument"+ | flag == "--out" = case rest of+ (v:rs) -> parseQuantileOpts rs (acc { qoOut = v })+ [] -> Left "--out requires a file path"+ | flag == "--report" = case rest of+ (v:rs) | not (null v) && head v /= '-' ->+ parseQuantileOpts rs (acc { qoReport = Just v })+ _ -> parseQuantileOpts rest (acc { qoReport = Just "quantile.html" })+ | otherwise = Left ("unexpected argument '" ++ flag ++ "'")++splitOnComma :: String -> [String]+splitOnComma s = case break (== ',') s of+ (a, ',' : rest) -> a : splitOnComma rest+ (a, _) -> [a]++doQuantile :: FilePath -> String -> String -> QuantileOpts -> LoadOpts -> IO ()+doQuantile file xColsStr yColStr opts lopts = do+ let xCols = map T.pack (words xColsStr)+ yCol = T.pack yColStr+ result <- loadXY lopts file xCols yCol+ case result of+ Left err -> hPutStrLn stderr err+ Right (df, xVecs, yVec) -> do+ let n = V.length yVec+ intercept = LA.konst 1 n+ xMat = LA.fromColumns+ (intercept : map (LA.fromList . V.toList) xVecs)+ yLA = LA.fromList (V.toList yVec)+ tau = qoTau opts+ fit = QR.fitQuantile tau xMat yLA+ beta = LA.toList (QR.qfBeta fit)+ printf "Loaded %d rows from %s\n" n file+ printf "Quantile: tau = %.3f (median: %s)\n" tau+ (if abs (tau - 0.5) < 1e-9 then ("yes" :: String) else "no")+ printf "MM-IRLS converged in %d iterations\n" (QR.qfIters fit)+ putStrLn ""+ putStrLn "Coefficients:"+ printf " %-30s = %9.4f\n" ("intercept" :: String) (head beta)+ mapM_ (\(i, c, b) ->+ printf " %-30s = %9.4f\n"+ ("β_" ++ show (i :: Int) ++ " (" ++ T.unpack c ++ ")") b)+ (zip3 [1..] xCols (tail beta))+ printf "Pinball loss V̂_τ: %.4f\n" (QR.qfPinball fit)+ printf "Pseudo R¹_τ: %.4f\n" (QR.qfR1 fit)++ -- 単変数なら scatter + fit (+ overlay multiple quantiles)+ case (xCols, xVecs) of+ ([xc1], [xVec]) -> do+ let xs = V.toList xVec+ ys = V.toList yVec+ grid = makeGrid xVec 100+ gridMat = LA.fromColumns+ [ LA.konst 1 100, LA.fromList grid ]+ gridY = LA.toList (QR.predictQuantile fit gridMat)+ sf = SmoothFit+ { sfX = grid+ , sfFit = gridY+ , sfLower = []+ , sfUpper = []+ , sfHasBand = False+ }+ _ = xs+ writeSmoothPlot (qoFormat opts) (qoOut opts)+ (T.pack ("Quantile τ=" ++ show tau)) df xc1 yCol sf+ putStrLn ("Plot: " ++ qoOut opts)+ openInBrowser (qoOut opts)++ -- HTML レポート+ case qoReport opts of+ Nothing -> return ()+ Just rpath -> do+ let coeffPairs = zip ("intercept" : xCols) beta+ modelLbl = "Quantile regression (τ=" <> T.pack (show tau) <> ")"+ formula = T.pack ("Q_τ(" ++ T.unpack yCol ++ "|x) = "+ ++ "β₀ + " ++ T.unpack (T.intercalate " + "+ [ "β" <> T.pack (show i)+ <> "·" <> c+ | (i, c) <- zip [(1::Int)..] xCols ]))+ cfg = RB.defaultReportConfig+ ("Quantile regression — τ=" <> T.pack (show tau)+ <> ", " <> yCol <> " ~ " <> T.intercalate " + " xCols)+ baseSections =+ [ RB.secDataOverview df xCols yCol+ , RB.secModelOverview modelLbl formula Nothing+ , RB.secCoefficients coeffPairs (Just ("Pseudo R¹_τ", QR.qfR1 fit))+ , RB.secKeyValue "Fit summary"+ [ ("τ", T.pack (printf "%.3f" tau))+ , ("Pinball loss V̂_τ", T.pack (printf "%.4f" (QR.qfPinball fit)))+ , ("Iterations", T.pack (show (QR.qfIters fit)))+ ]+ , RB.secFitScatter xc1 yCol xs ys (Just (RB.SmoothCurve grid gridY [] []))+ , RB.secResiduals (LA.toList (QR.qfYHat fit))+ (LA.toList (QR.qfResid fit))+ ]+ -- overlay multi quantile chart+ multiSec = case qoTaus opts of+ [] -> []+ taus ->+ let curves = [ ( T.pack ("τ=" ++ show t)+ , LA.toList (QR.predictQuantile+ (QR.fitQuantile t xMat yLA)+ gridMat))+ | t <- taus ]+ spec = multiQuantileSpec xc1 yCol xs ys grid curves+ in [RB.secVega "Multiple quantile fits" spec]+ RB.renderReport rpath cfg (baseSections ++ multiSec)+ putStrLn ("Report: " ++ rpath)+ openInBrowser rpath+ _ -> putStrLn "(scatter plot skipped for multiple x columns)"++-- 複数分位線を 1 枚の Vega-Lite spec で描く+multiQuantileSpec :: T.Text -> T.Text -> [Double] -> [Double] -> [Double]+ -> [(T.Text, [Double])] -> VegaLite+multiQuantileSpec xc yc xs ys grid curves =+ VL.toVegaLite+ [ VL.layer+ [ VL.asSpec+ [ VL.dataFromColumns []+ . VL.dataColumn xc (VL.Numbers xs)+ . VL.dataColumn yc (VL.Numbers ys)+ $ []+ , VL.mark VL.Point+ [VL.MOpacity 0.5, VL.MSize 40, VL.MColor "#888888"]+ , VL.encoding+ . VL.position VL.X+ [VL.PName xc, VL.PmType VL.Quantitative,+ VL.PAxis [VL.AxTitle xc]]+ . VL.position VL.Y+ [VL.PName yc, VL.PmType VL.Quantitative,+ VL.PAxis [VL.AxTitle yc]]+ $ []+ ]+ , VL.asSpec (multiLineLayer xc yc grid curves)+ ]+ , VL.width 640+ , VL.height 320+ ]++multiLineLayer :: T.Text -> T.Text -> [Double] -> [(T.Text, [Double])]+ -> [(VLProperty, VLSpec)]+multiLineLayer xc yc grid curves =+ let rowsX = concat [ replicate (length grid) lbl | (lbl, _) <- curves ]+ rowsXs = concat [ grid | _ <- curves ]+ rowsYs = concat [ ys' | (_, ys') <- curves ]+ in [ VL.dataFromColumns []+ . VL.dataColumn "tau" (VL.Strings rowsX)+ . VL.dataColumn xc (VL.Numbers rowsXs)+ . VL.dataColumn yc (VL.Numbers rowsYs)+ $ []+ , VL.mark VL.Line [VL.MStrokeWidth 2.2]+ , VL.encoding+ . VL.position VL.X [VL.PName xc, VL.PmType VL.Quantitative]+ . VL.position VL.Y [VL.PName yc, VL.PmType VL.Quantitative]+ . VL.color [VL.MName "tau", VL.MmType VL.Nominal,+ VL.MScale [VL.SScheme "tableau10" []]]+ $ []+ ]++-- ---------------------------------------------------------------------------+-- gam subcommand+-- ---------------------------------------------------------------------------++gamUsage :: String+gamUsage = unlines+ [ "Usage: hanalyze gam <file> <xcols> <ycol> [options]"+ , ""+ , " <xcols> x column names; quote multiple: \"x1 x2 x3\""+ , " <ycol> y column name"+ , ""+ , "Options:"+ , " --knots N per-feature internal knot count (default: 5)"+ , " --degree D B-spline degree (default: 3 = cubic)"+ , " --lambda L Ridge regularization on spline coefficients (default: 0.01)"+ , " --report [FILE] build composite HTML report with per-feature partials"+ , ""+ , "Example:"+ , " hanalyze gam data.csv \"x1 x2 x3\" y --knots 8 --lambda 0.05 --report"+ ]++data GAMOpts = GAMOpts+ { goKnots :: Int+ , goDegree :: Int+ , goLambda :: Double+ , goReport :: Maybe FilePath+ }++defaultGAMOpts :: GAMOpts+defaultGAMOpts = GAMOpts 5 3 0.01 Nothing++runGAMCmd :: [String] -> IO ()+runGAMCmd args0 =+ let (lopts, args) = parseLoadOpts args0+ in case args of+ (file : xColsStr : yColStr : rest) ->+ case parseGAMOpts rest defaultGAMOpts of+ Left err -> hPutStrLn stderr ("gam: " ++ err)+ Right opts -> doGAM file xColsStr yColStr opts lopts+ _ -> putStrLn gamUsage++parseGAMOpts :: [String] -> GAMOpts -> Either String GAMOpts+parseGAMOpts [] acc = Right acc+parseGAMOpts (flag:rest) acc+ | flag == "--knots" = case rest of+ (v:rs) -> case reads v :: [(Int,String)] of+ [(d,"")] -> parseGAMOpts rs (acc { goKnots = d })+ _ -> Left ("invalid --knots '" ++ v ++ "'")+ [] -> Left "--knots requires a value"+ | flag == "--degree" = case rest of+ (v:rs) -> case reads v :: [(Int,String)] of+ [(d,"")] -> parseGAMOpts rs (acc { goDegree = d })+ _ -> Left ("invalid --degree '" ++ v ++ "'")+ [] -> Left "--degree requires a value"+ | flag == "--lambda" = case rest of+ (v:rs) -> case reads v :: [(Double,String)] of+ [(d,"")] -> parseGAMOpts rs (acc { goLambda = d })+ _ -> Left ("invalid --lambda '" ++ v ++ "'")+ [] -> Left "--lambda requires a value"+ | flag == "--report" = case rest of+ (v:rs) | not (null v) && head v /= '-' ->+ parseGAMOpts rs (acc { goReport = Just v })+ _ -> parseGAMOpts rest (acc { goReport = Just "gam.html" })+ | otherwise = Left ("unexpected argument '" ++ flag ++ "'")++doGAM :: FilePath -> String -> String -> GAMOpts -> LoadOpts -> IO ()+doGAM file xColsStr yColStr opts lopts = do+ let xCols = map T.pack (words xColsStr)+ yCol = T.pack yColStr+ result <- loadXY lopts file xCols yCol+ case result of+ Left err -> hPutStrLn stderr err+ Right (df, xVecs, yVec) -> do+ let fit = GAM.fitGAM (goDegree opts) (goKnots opts) (goLambda opts)+ xVecs yVec+ n = V.length yVec+ ys = V.toList yVec+ yhat = LA.toList (GAM.gamYHat fit)+ resid = LA.toList (GAM.gamResid fit)+ printf "Loaded %d rows from %s\n" n file+ printf "GAM: degree=%d, knots=%d/feature, lambda=%g\n"+ (goDegree opts) (goKnots opts) (goLambda opts)+ printf "Features: %d (%s)\n" (length xCols)+ (T.unpack (T.intercalate ", " xCols))+ printf "Intercept: %.4f\n" (GAM.gamIntercept fit)+ printf "R²: %.4f\n" (GAM.gamR2 fit)+ let rmseVal = sqrt (sum [ r ^ (2 :: Int) | r <- resid ]+ / fromIntegral n)+ printf "RMSE (in-sample): %.4f\n" rmseVal++ case goReport opts of+ Nothing -> return ()+ Just rpath -> do+ let modelLbl = "Generalized Additive Model"+ formula = yCol <> " = β₀ + " <> T.intercalate " + "+ [ "s(" <> c <> ")" | c <- xCols ]+ cfg = RB.defaultReportConfig+ ("GAM — " <> yCol <> " ~ s("+ <> T.intercalate ") + s(" xCols <> ")")+ partialSecs =+ [ RB.secVega ("Partial effect: s(" <> c <> ")")+ (gamPartialSpec c xVec fit j)+ | (j, c, xVec) <- zip3 [0..] xCols xVecs ]+ sections =+ [ RB.secDataOverview df xCols yCol+ , RB.secModelOverview modelLbl formula Nothing+ , RB.secKeyValue "Fit summary"+ [ ("Degree", T.pack (show (goDegree opts)))+ , ("Knots", T.pack (show (goKnots opts)))+ , ("Lambda", T.pack (printf "%g" (goLambda opts)))+ , ("Intercept",T.pack (printf "%.4f"+ (GAM.gamIntercept fit)))+ , ("R²", T.pack (printf "%.4f" (GAM.gamR2 fit)))+ , ("RMSE", T.pack (printf "%.4f" rmseVal))+ ]+ ] ++ partialSecs +++ [ RB.secResiduals yhat resid ]+ _ = ys+ RB.renderReport rpath cfg sections+ putStrLn ("Report: " ++ rpath)+ openInBrowser rpath++-- ---------------------------------------------------------------------------+-- rf subcommand+-- ---------------------------------------------------------------------------++rfUsage :: String+rfUsage = unlines+ [ "Usage: hanalyze rf <file> <xcols> <ycol> [options]"+ , ""+ , "Options:"+ , " --trees N number of trees (default: 100)"+ , " --max-depth D maximum tree depth (default: 12)"+ , " --min-samples N minimum samples per leaf (default: 3)"+ , " --mtry M features per split (default: max(1, d/3))"+ , " --report [FILE] build composite HTML report (with feature importance)"+ , ""+ , "Example:"+ , " hanalyze rf data.csv \"x1 x2 x3\" y --trees 200 --report"+ ]++data RFOpts = RFOpts+ { roTrees :: Int+ , roMaxDepth :: Int+ , roMinSamples :: Int+ , roMtry :: Maybe Int+ , roReport_ :: Maybe FilePath+ }++defaultRFOpts :: RFOpts+defaultRFOpts = RFOpts 100 12 3 Nothing Nothing++runRFCmd :: [String] -> IO ()+runRFCmd args0 =+ let (lopts, args) = parseLoadOpts args0+ in case args of+ (file : xColsStr : yColStr : rest) ->+ case parseRFOpts rest defaultRFOpts of+ Left err -> hPutStrLn stderr ("rf: " ++ err)+ Right opts -> doRF file xColsStr yColStr opts lopts+ _ -> putStrLn rfUsage++parseRFOpts :: [String] -> RFOpts -> Either String RFOpts+parseRFOpts [] acc = Right acc+parseRFOpts (flag:rest) acc+ | flag == "--trees" = case rest of+ (v:rs) -> case reads v :: [(Int,String)] of+ [(d,"")] -> parseRFOpts rs (acc { roTrees = d })+ _ -> Left ("invalid --trees '" ++ v ++ "'")+ [] -> Left "--trees requires a value"+ | flag == "--max-depth" = case rest of+ (v:rs) -> case reads v :: [(Int,String)] of+ [(d,"")] -> parseRFOpts rs (acc { roMaxDepth = d })+ _ -> Left ("invalid --max-depth '" ++ v ++ "'")+ [] -> Left "--max-depth requires a value"+ | flag == "--min-samples" = case rest of+ (v:rs) -> case reads v :: [(Int,String)] of+ [(d,"")] -> parseRFOpts rs (acc { roMinSamples = d })+ _ -> Left ("invalid --min-samples '" ++ v ++ "'")+ [] -> Left "--min-samples requires a value"+ | flag == "--mtry" = case rest of+ (v:rs) -> case reads v :: [(Int,String)] of+ [(d,"")] -> parseRFOpts rs (acc { roMtry = Just d })+ _ -> Left ("invalid --mtry '" ++ v ++ "'")+ [] -> Left "--mtry requires a value"+ | flag == "--report" = case rest of+ (v:rs) | not (null v) && head v /= '-' ->+ parseRFOpts rs (acc { roReport_ = Just v })+ _ -> parseRFOpts rest (acc { roReport_ = Just "rf.html" })+ | otherwise = Left ("unexpected argument '" ++ flag ++ "'")++doRF :: FilePath -> String -> String -> RFOpts -> LoadOpts -> IO ()+doRF file xColsStr yColStr opts lopts = do+ let xCols = map T.pack (words xColsStr)+ yCol = T.pack yColStr+ result <- loadXY lopts file xCols yCol+ case result of+ Left err -> hPutStrLn stderr err+ Right (df, xVecs, yVec) -> do+ let n = V.length yVec+ rows = [ [ xv V.! i | xv <- xVecs ] | i <- [0 .. n - 1] ]+ ys = V.toList yVec+ cfg = RF.defaultRandomForest+ { RF.rfTrees = roTrees opts+ , RF.rfMaxDepth = roMaxDepth opts+ , RF.rfMinSamples = roMinSamples opts+ , RF.rfMtry = roMtry opts+ }+ gen <- createSystemRandom+ forest <- RF.fitRF cfg rows ys gen+ let yhat = map (RF.predictRF forest) rows+ resid = zipWith (-) ys yhat+ yMean = sum ys / fromIntegral n+ tss = sum [ (y - yMean) ^ (2 :: Int) | y <- ys ]+ rss = sum [ r ^ (2 :: Int) | r <- resid ]+ r2 = if tss < 1e-12 then 0 else 1 - rss / tss+ rmseVal = sqrt (rss / fromIntegral n)+ imp = V.toList (RF.featureImportance forest)+ impPairs = zip xCols imp+ printf "Loaded %d rows from %s\n" n file+ printf "RandomForest: trees=%d, max-depth=%d, min-samples=%d\n"+ (roTrees opts) (roMaxDepth opts) (roMinSamples opts)+ printf "R²: %.4f\n" r2+ printf "RMSE (in-sample): %.4f\n" rmseVal+ putStrLn ""+ putStrLn "Feature importance (split-count fraction):"+ mapM_ (\(c, v) -> printf " %-20s = %.4f\n" (T.unpack c) v) impPairs++ case roReport_ opts of+ Nothing -> return ()+ Just rpath -> do+ let modelLbl = "Random Forest regression"+ formula = yCol <> " ~ ensemble of " <> T.pack (show (roTrees opts))+ <> " CART trees over (" <> T.intercalate ", " xCols <> ")"+ cfg' = RB.defaultReportConfig+ ("Random Forest — " <> yCol <> " ~ "+ <> T.intercalate " + " xCols)+ sections =+ [ RB.secDataOverview df xCols yCol+ , RB.secModelOverview modelLbl formula Nothing+ , RB.secKeyValue "Fit summary"+ [ ("Trees", T.pack (show (roTrees opts)))+ , ("Max depth", T.pack (show (roMaxDepth opts)))+ , ("Min samples", T.pack (show (roMinSamples opts)))+ , ("R²", T.pack (printf "%.4f" r2))+ , ("RMSE", T.pack (printf "%.4f" rmseVal))+ ]+ , RB.secBarChart "Feature importance"+ [ (c, v) | (c, v) <- impPairs ]+ , RB.secResiduals yhat resid+ ]+ RB.renderReport rpath cfg' sections+ putStrLn ("Report: " ++ rpath)+ openInBrowser rpath++-- 1 特徴の partial effect s_j(x_j) を Vega-Lite 散布+曲線で+gamPartialSpec :: T.Text -> V.Vector Double -> GAM.GAMFit -> Int -> VegaLite+gamPartialSpec col xVec fit j =+ let xs = V.toList xVec+ lo = V.minimum xVec+ hi = V.maximum xVec+ grid = [ lo + fromIntegral i * (hi - lo) / 99 | i <- [0..99::Int]]+ gridV = V.fromList grid+ sj = V.toList (GAM.predictGAMComponent fit j gridV)+ -- partial residuals: resid + s_j(x_i) (説明用にプロット)+ partialAtData = V.toList (GAM.predictGAMComponent fit j xVec)+ residList = LA.toList (GAM.gamResid fit)+ partials = zipWith (+) residList partialAtData+ in VL.toVegaLite+ [ VL.layer+ [ VL.asSpec+ [ VL.dataFromColumns []+ . VL.dataColumn col (VL.Numbers xs)+ . VL.dataColumn "partial" (VL.Numbers partials)+ $ []+ , VL.mark VL.Point+ [VL.MOpacity 0.5, VL.MSize 40, VL.MColor "#888888"]+ , VL.encoding+ . VL.position VL.X+ [VL.PName col, VL.PmType VL.Quantitative,+ VL.PAxis [VL.AxTitle col]]+ . VL.position VL.Y+ [VL.PName "partial", VL.PmType VL.Quantitative,+ VL.PAxis [VL.AxTitle "Partial residual"]]+ $ []+ ]+ , VL.asSpec+ [ VL.dataFromColumns []+ . VL.dataColumn col (VL.Numbers grid)+ . VL.dataColumn "s_j" (VL.Numbers sj)+ $ []+ , VL.mark VL.Line+ [VL.MStrokeWidth 2.5, VL.MColor "#DD5566"]+ , VL.encoding+ . VL.position VL.X+ [VL.PName col, VL.PmType VL.Quantitative]+ . VL.position VL.Y+ [VL.PName "s_j", VL.PmType VL.Quantitative]+ $ []+ ]+ ]+ , VL.width 500+ , VL.height 240+ ]++-- ---------------------------------------------------------------------------+-- multireg subcommand (多出力回帰: wide CSV → 対話的予測曲線)+-- ---------------------------------------------------------------------------++multiRegUsage :: String+multiRegUsage = unlines+ [ "Usage: hanalyze multireg <file> <xcol> <yspec> [options]"+ , ""+ , "wide-form CSV (1 行 = 入力 1 値、複数列 = q 個の出力) を読み込み、"+ , "1 入力 → q 出力の多出力回帰を実行。dose スライダで対話的に予測曲線を更新。"+ , ""+ , "<yspec>: カンマ区切り列名 (例 'y_z001,y_z002,...') または prefix*"+ , " (例 'y_z*' で y_z で始まる全列)"+ , ""+ , "Options:"+ , " --method M linear | kernel-rbf (default: linear)"+ , " --bandwidth H kernel-rbf の bandwidth (default: auto via LOOCV)"+ , " --lambda L kernel-rbf の Ridge λ (default: auto via LOOCV)"+ , " --auto-hp kernel-rbf で h, λ を LOOCV 解析解で自動決定 (default: ON)"+ , " --report FILE 対話的 HTML レポート出力先 (default: multireg.html)"+ , " --xaxis LABEL 出力グリッドの x 軸ラベル (default: 'index')"+ , ""+ , "前提: y 列が共通の z grid を表す場合、列名末尾の数値で z 座標を内挿。"+ , " 例: y_z001..y_z100 のとき z = 0..99 を等間隔展開 (--xaxis-min/max で上書き)."+ , ""+ , "Options (出力 grid):"+ , " --xaxis-min V 出力 grid の最小値 (default: 1)"+ , " --xaxis-max V 出力 grid の最大値 (default: q)"+ , ""+ , "Examples:"+ , " hanalyze multireg data/io/potential_wide.csv dose 'y_z*' \\"+ , " --method kernel-rbf --report trash/pot.html \\"+ , " --xaxis 'z [nm]' --xaxis-min 0 --xaxis-max 200"+ ]++data MROpts = MROpts+ { mroMethod :: String -- "linear" | "kernel-rbf"+ , mroH :: Maybe Double+ , mroLambda :: Maybe Double+ , mroAutoHP :: Bool+ , mroReport :: FilePath+ , mroXAxis :: String+ , mroXAxisMin :: Maybe Double+ , mroXAxisMax :: Maybe Double+ } deriving Show++defaultMROpts :: MROpts+defaultMROpts = MROpts "linear" Nothing Nothing True "multireg.html" "index" Nothing Nothing++parseMROpts :: [String] -> (MROpts, [String])+parseMROpts = go defaultMROpts []+ where+ go o acc [] = (o, reverse acc)+ go o acc ("--method":m:rest) = go o { mroMethod = m } acc rest+ go o acc ("--bandwidth":v:rest) = go o { mroH = Just (read v) } acc rest+ go o acc ("--lambda":v:rest) = go o { mroLambda = Just (read v) } acc rest+ go o acc ("--auto-hp":rest) = go o { mroAutoHP = True } acc rest+ go o acc ("--no-auto-hp":rest) = go o { mroAutoHP = False } acc rest+ go o acc ("--report":p:rest) = go o { mroReport = p } acc rest+ go o acc ("--xaxis":s:rest) = go o { mroXAxis = s } acc rest+ go o acc ("--xaxis-min":v:rest) = go o { mroXAxisMin = Just (read v) } acc rest+ go o acc ("--xaxis-max":v:rest) = go o { mroXAxisMax = Just (read v) } acc rest+ go o acc (x:rest) = go o (x:acc) rest++runMultiRegCmd :: [String] -> IO ()+runMultiRegCmd args0 = do+ let (lopts, args1) = parseLoadOpts args0+ (opts, args2) = parseMROpts args1+ case args2 of+ (file:xCol:ySpec:_) -> do+ result <- loadAutoSafeWith lopts file+ case result of+ Left err -> hPutStrLn stderr ("Parse error: " ++ err)+ Right (df, lg) -> do+ Log.printLogReport lg+ let allCols = map T.unpack (DX.columnNames df)+ yCols = resolveYSpec ySpec allCols+ xColT = T.pack xCol+ yColTs = map T.pack yCols+ if null yCols+ then hPutStrLn stderr ("multireg: yspec '" ++ ySpec+ ++ "' に該当する列がありません")+ else case getDoubleVec xColT df of+ Nothing -> hPutStrLn stderr ("multireg: 入力列 '" ++ xCol+ ++ "' が見つかりません")+ Just xV -> do+ let n = V.length xV+ yMs = [ getDoubleVec c df | c <- yColTs ]+ if any null (map mtoMaybe yMs)+ then hPutStrLn stderr "multireg: y 列の取得失敗"+ else do+ let yVecs = [v | Just v <- yMs]+ q = length yVecs+ xMat1 = LA.fromLists [[1.0, xV V.! i]+ | i <- [0 .. n - 1]]+ ys = LA.fromLists+ [ [ (yVecs !! j) V.! i+ | j <- [0 .. q - 1] ]+ | i <- [0 .. n - 1] ]+ xObsL = V.toList xV+ yObsL = [ [ (yVecs !! j) V.! i+ | j <- [0 .. q - 1] ]+ | i <- [0 .. n - 1] ]+ outGrid =+ let lo = maybe 1.0 id (mroXAxisMin opts)+ hi = maybe (fromIntegral q) id (mroXAxisMax opts)+ step = if q < 2 then 0 else (hi - lo) / fromIntegral (q - 1)+ in [ lo + step * fromIntegral i | i <- [0 .. q - 1] ]+ xMin = minimum xObsL - (maximum xObsL - minimum xObsL) * 0.2+ xMax = maximum xObsL + (maximum xObsL - minimum xObsL) * 0.2+ xMid = 0.5 * (xMin + xMax)+ putStrLn $ "Loaded " ++ show n ++ " rows × " ++ show q+ ++ " outputs; method=" ++ mroMethod opts+ sections <- case mroMethod opts of+ "linear" -> do+ let mf = MLM.fitMultiLM xMat1 ys+ betaB = Core.coefficients (MLM.mfFit mf)+ ints = LA.toList (betaB LA.! 0)+ slps = LA.toList (betaB LA.! 1)+ res = Core.residuals (MLM.mfFit mf)+ rmse = sqrt (LA.sumElements (res*res)+ / fromIntegral (n * q))+ r2v = Core.rSquared (MLM.mfFit mf)+ r2mean = LA.sumElements r2v / fromIntegral q+ imo = RB.mkInteractiveMOLinear+ (T.pack xCol)+ "y" (T.pack (mroXAxis opts))+ outGrid xObsL yObsL+ ints slps (xMin, xMid, xMax)+ printf " RMSE = %.4f, R^2 mean = %.4f\n" rmse r2mean+ return+ [ RB.secModelOverview "Multi-output Linear (B = (X'X)^-1 X'Y)"+ "$\\hat{Y} = X B$" Nothing+ , RB.secStatRow+ [ ("N", T.pack (show n))+ , ("q", T.pack (show q))+ , ("RMSE", T.pack (printf "%.4f" rmse))+ , ("R^2 mean", T.pack (printf "%.4f" r2mean))+ ]+ , RB.secInteractiveMultiOut "予測曲線 (スライダ)" imo+ ]+ "kernel-rbf" -> do+ let hs = case mroH opts of+ Just h | not (mroAutoHP opts) -> [h]+ _ -> Kern.defaultHGrid xV+ lams = case mroLambda opts of+ Just l | not (mroAutoHP opts) -> [l]+ _ -> Kern.defaultLamGrid+ (fit, bestH, bestL, looMSE) =+ Kern.autoTuneKernelRidgeMulti+ Kern.Gaussian xV ys hs lams+ yhat = Kern.fittedKernelRidgeMulti fit+ r2v = Kern.r2Multi ys yhat+ res = ys - yhat+ rmse = sqrt (LA.sumElements (res*res)+ / fromIntegral (n * q))+ r2mean = V.sum r2v / fromIntegral q+ alpha2 = [ LA.toList (LA.flatten (Kern.krmAlpha fit LA.? [i]))+ | i <- [0 .. n - 1] ]+ imo = RB.mkInteractiveMOKernelRBF+ (T.pack xCol) "y"+ (T.pack (mroXAxis opts))+ outGrid xObsL yObsL+ xObsL alpha2 bestH+ (xMin, xMid, xMax)+ printf " best h=%.3g λ=%.3g LOO MSE=%.3g RMSE=%.4f R^2 mean=%.4f\n"+ bestH bestL looMSE rmse r2mean+ return+ [ RB.secModelOverview "Multi-output Kernel Ridge (RBF)"+ "$\\hat{y}_j(x)=\\sum_i K_h(x,x_i)\\,\\alpha_{ij}$" Nothing+ , RB.secStatRow+ [ ("N", T.pack (show n))+ , ("q", T.pack (show q))+ , ("h", T.pack (printf "%.3g" bestH))+ , ("λ", T.pack (printf "%.3g" bestL))+ , ("LOO MSE", T.pack (printf "%.3g" looMSE))+ , ("RMSE", T.pack (printf "%.4f" rmse))+ , ("R^2 mean", T.pack (printf "%.4f" r2mean))+ ]+ , RB.secInteractiveMultiOut "予測曲線 (スライダ)" imo+ ]+ m -> do+ hPutStrLn stderr ("multireg: unknown method '" ++ m+ ++ "' (use linear|kernel-rbf)")+ return []+ if null sections+ then return ()+ else do+ let cfg = RB.defaultReportConfig+ (T.pack ("multireg: " ++ file))+ RB.renderReport (mroReport opts) cfg sections+ putStrLn ("Wrote " ++ mroReport opts)+ _ -> hPutStrLn stderr multiRegUsage+ where+ mtoMaybe Nothing = []+ mtoMaybe (Just _) = ["x"]+ -- "y_z*" → all columns starting with "y_z"+ -- "a,b,c" → ["a","b","c"]+ resolveYSpec spec allCols+ | last' spec == Just '*' =+ let pre = init spec+ in [ c | c <- allCols, take (length pre) c == pre, c /= xCol0 spec ]+ | otherwise = wordsBy (== ',') spec+ last' [] = Nothing+ last' s = Just (last s)+ -- xCol0 is irrelevant for filtering but ensure no accidental match+ xCol0 _ = ""+
+ hanalyze-cli.cabal view
@@ -0,0 +1,53 @@+cabal-version: 3.0+name: hanalyze-cli+version: 0.2.0.1+synopsis: hanalyze command-line interface for the hanalyze toolkit+description:+ The @hanalyze@ command-line executable of the hanalyze toolkit.+ Its 15 subcommands cover regression (@regress@ for LM / GLM / GLMM / GP /+ HBM, plus @ridge@, @kernel@, @spline@, @quantile@, @gam@, @rf@ and+ @multireg@), data inspection and plotting (@info@, @hist@), design of+ experiments (@doe@, @taguchi@) and reshaping (@clean@, @melt@, @regrid@),+ with HTML / PNG / SVG output and optional HTML analysis reports.+ .+ It was split out of the hanalyze package so that library work does+ not trigger the CLI compile+link by default, and so that its dependency on+ the umbrella public API alone (no internal module access) is guaranteed+ structurally. See README.md for the subcommand map and examples.+license: BSD-3-Clause+author: Toshiaki Honda+maintainer: frenzieddoll@gmail.com+copyright: 2026 Aelysce Project (Toshiaki Honda)+category: Math, Statistics, Numeric, Machine Learning+build-type: Simple+tested-with: GHC == 9.6.7+extra-source-files:+ README.md+ README.ja.md++common warnings+ ghc-options: -Wall -Wcompat -Widentities -Wredundant-constraints++common opt+ ghc-options: -O2 -funbox-strict-fields++-- exe 名は分割前と同じ hanalyze を維持 (cabal run hanalyze -- ... 互換)+executable hanalyze+ import: warnings, opt+ main-is: Main.hs+ hs-source-dirs: app+ default-language: GHC2021+ build-depends:+ base >= 4.14 && < 5+ , hanalyze == 0.2.0.1+ , text >= 1.2 && < 2.2+ , vector >= 0.12 && < 0.14+ , hmatrix >= 0.20 && < 0.22+ , mwc-random >= 0.15 && < 0.16+ , containers >= 0.6 && < 0.8+ , filepath >= 1.4 && < 1.6+ , hvega >= 0.12 && < 0.13+ , dataframe-core ^>= 1.1+ , dataframe-operations >= 1.1.1 && < 1.2+ , dataframe-csv ^>= 1.0.2+ , time >= 1.11 && < 1.13