using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Windows.Forms;

namespace vcsMRA
{
    public partial class Form1 : Form
    {
        public Form1()
        {
            InitializeComponent();
        }
        //--------------------------------------------------------------------------
        private void toolStripButton1_Click(object sender, EventArgs e)
        {
            //重回帰分析（multiple regression analysis）
            //入出力用変数
            string fnameR = "";
            string fnameW = "";
            System.IO.StreamReader sr;
            System.IO.StreamWriter sw;
            string dat="";
            string[] sbuf;
            char delim = ',';

            string scom;   //コメント
            int ndata;     //データ組数−１
            int mcol;      //説明変数の数−１
            int mm;        //説明変数の数−１＋１
            double RR;     //重相関係数
            double YEmean; //回帰推定値の平均
            double YDmean; //目的変数入力値の平均
            double std;    //残差の標準偏差

            int i, j, k;
            double y1, y2, y3;

            listBox1.Items.Clear();

            //入出力ファイル指定
            openFileDialog1.InitialDirectory=System.IO.Directory.GetCurrentDirectory();
            if(openFileDialog1.ShowDialog()==DialogResult.OK){fnameR=openFileDialog1.FileName;}
            if(saveFileDialog1.ShowDialog()==DialogResult.OK){fnameW=saveFileDialog1.FileName;}

	        //データ読み込み
    	    sr = new System.IO.StreamReader(fnameR, System.Text.Encoding.Default);
        	dat = sr.ReadLine() ; sbuf = dat.Split(delim) ; scom = sbuf[0];
        	dat = sr.ReadLine() ; sbuf = dat.Split(delim) ; mcol = int.Parse(sbuf[0]) - 1 ; ndata = int.Parse(sbuf[1]) - 1;
        	mm = mcol + 1;
        	double[,] XD=new double[ndata+1, mm+1]; //説明変数入力値用配列
        	double[,] XT=new double[mm+1, ndata+1]; //説明変数入力値用配列の転置行列
        	double[]  YD=new double[ndata+1];       //目的変数入力値
        	double[,] a=new double[mm+1, mm + 2];   //連立方程式用配列
        	double[]  BB=new double[mm+1];          //偏回帰係数
        	double[]  YE=new double[ndata+1];       //回帰推定値
        	double[]  RS=new double[ndata+1];       //目的変数入力値と回帰推定値の差分
        	i = 0;
           	while(!sr.EndOfStream){
	            dat = sr.ReadLine() ; sbuf = dat.Split(delim);
            	XD[i, 0] = 1.0;
            	for(j = 0;j<=mcol;j++){
                	XD[i, j + 1] = double.Parse(sbuf[j]);
            	}
            	YD[i] = double.Parse(sbuf[mcol + 1]);
            	i = i + 1;
            }
        	sr.Close();

            //正規方程式作成
            for(i = 0;i<=ndata;i++){
                for(j = 0;j<=mm;j++){
                    XT[j, i] = XD[i, j];
            }}
            for(i = 0;i<=mm;i++){
                for(j = 0;j<=mm;j++){
                    a[i, j] = 0.0;
                    for(k = 0;k<=ndata;k++){
                        a[i, j] = a[i, j] + XT[i, k] * XD[k, j];
                    }
		    }}
            for(i = 0;i<=mm;i++){
                a[i, mm + 1] = 0.0;
                for(k = 0;k<=ndata;k++){
                    a[i, mm + 1] = a[i, mm + 1] + XT[i, k] * YD[k];
                }
            }

            //連立一次方程式の解
            MATGJ(mm, a);
            for(i = 0;i<=mm;i++){
                BB[i] = a[i, mm + 1];
            }

            //重相関係数算定
            YDmean = 0.0 ; YEmean = 0.0;
            for(i = 0;i<=ndata;i++){
                YE[i] = 0.0;
                for(j = 0;j<=mm;j++){
                    YE[i] = YE[i] + XD[i, j] * BB[j];
                }
                YDmean = YDmean + YD[i];
                YEmean = YEmean + YE[i];
            }
            YDmean = YDmean / (double)(ndata + 1);
            YEmean = YEmean / (double)(ndata + 1);
            y1 = 0.0 ; y2 = 0.0 ; y3 = 0.0;
            for(i = 0;i<=ndata;i++){
                y1 = y1 + (YD[i] - YDmean) * (YE[i] - YEmean);
                y2 = y2 + (YD[i] - YDmean) * (YD[i] - YDmean);
                y3 = y3 + (YE[i] - YEmean) * (YE[i] - YEmean);
            }
            RR = y1 / Math.Sqrt(y2 * y3);
            //残差の標準偏差
            y1 = 0.0;
            for(i = 0;i<=ndata;i++){
                RS[i] = YD[i] - YE[i];
                y1 = y1 + RS[i];
            }
            y1 = y1 / (double)(ndata + 1); //残差の平均
            y2 = 0.0;
            for(i = 0;i<=ndata;i++){
                y2 = y2 + (RS[i] - y1)*(RS[i] - y1);
            }
            std = Math.Sqrt(y2 / (double)(ndata)); //残差の標準偏差

            //分析結果の画面（ListBox）への書き出し
            listBox1.Items.Add(scom);
            listBox1.Items.Add("y=B0+B1*x1+B2*x2+・・・Bm*xm");
            for(i = 0;i<=mm;i++){
                listBox1.Items.Add("B" + i.ToString("0") + "=" + BB[i].ToString("E"));
            }
            listBox1.Items.Add("重相関係数：R=" + RR.ToString("0.000"));
            listBox1.Items.Add("データ数：n=" + (ndata + 1).ToString("0"));
            listBox1.Items.Add("残差平均：" + y1.ToString("E"));
            listBox1.Items.Add("残差標準偏差：" + std.ToString("E"));

            //分析結果のファイルへの書き出し
            sw = new System.IO.StreamWriter(fnameW, false, System.Text.Encoding.Default);
            sw.WriteLine(scom);
            sw.WriteLine("y=B0+B1*x1+B2*x2+・・・Bm*xm");
            sw.WriteLine("偏回帰係数");
            for(i = 0;i<=mm;i++){
                sw.WriteLine("B" + i.ToString("0") + "=," + BB[i].ToString("E"));
            }
            sw.WriteLine("重相関係数：R=," + RR.ToString("0.000"));
            sw.WriteLine("データ数：n=," + (ndata + 1).ToString("0"));
            sw.WriteLine("残差平均：," + y1.ToString("E"));
            sw.WriteLine("残差標準偏差：," + std.ToString("E"));
            sw.WriteLine("目的変数入力値,回帰推定値");
            for(i = 0;i<=ndata;i++){
                sw.WriteLine(YD[i].ToString("E") + "," + YE[i].ToString("E"));
            }
            sw.Close();

        }
        //--------------------------------------------------------------------------
        private void MATGJ(int n, double[,] a)
        {
            //Gauss-Jordan法による連立一次方程式の解法
            int i, j, k, s;
            double p, d, max, dumy;

            for (k = 0; k <= n; k++)
            {
                max = 0.0;
                s = k;
                for (j = k; j <= n; j++)
                {
                    if (Math.Abs(a[j, k]) > max)
                    {
                        max = Math.Abs(a[j, k]);
                        s = j;
                    }
                }
                for (j = 0; j <= n + 1; j++)
                {
                    dumy = a[k, j];
                    a[k, j] = a[s, j];
                    a[s, j] = dumy;
                }
                p = a[k, k];
                for (j = k; j <= n + 1; j++)
                {
                    a[k, j] = a[k, j] / p;
                }
                for (i = 0; i <= n; i++)
                {
                    if (i != k)
                    {
                        d = a[i, k];
                        for (j = k; j <= n + 1; j++)
                        {
                            a[i, j] = a[i, j] - d * a[k, j];
                        }
                    }
                }
            }
        }
        //--------------------------------------------------------------------------

    }
}