-
Notifications
You must be signed in to change notification settings - Fork 7
/
矩阵运算类.cs
346 lines (315 loc) · 9.94 KB
/
矩阵运算类.cs
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
class CMat
{
//构建零矩阵
public double[,] MatZero(int m, int n)
{
double[,] result = new double[m, n];
int i, j;
for (i = 0; i <= m - 1; i++)
{
for (j = 0; j <= n - 1; j++)
{
result[i, j] = 0.0;
}
}
return result;
}
public double[] MatZero(int m)
{
double[] result = new double[m];
int i;
for (i = 0; i <= m - 1; i++)
{
result[i] = 0.0;
}
return result;
}
//乘-1
public double[] MatNegat(double[] a)
{
int m = a.GetLength(0);
double[] result = new double[m];
int i;
for (i = 0; i <= m - 1; i++)
{
result[i] = (-1) * a[i];
}
return result;
}
//求矩阵转置
public double[,] MatTrans(double[,] a)
{
double[,] result = null;
int m = a.GetLength(0);
int n = a.GetLength(1);
double[,] b = new double[n, m];
int i, j;
for (i = 0; i < m; i++)
{
for (j = 0; j < n; j++)
{
b[j, i] = a[i, j];
}
}
result = b;
return result;
}
//求矩阵相乘
public double[,] MatMulti(double[,] a, double[,] b)
{
double[,] result = null;
int m1 = a.GetLength(0);
int n1 = a.GetLength(1);
int m2 = b.GetLength(0);
int n2 = b.GetLength(1);
if (n1 != m2) return null;
double[,] c = new double[m1, n2];
int i, j, k;
for (i = 0; i < m1; i++)
{
for (j = 0; j < n2; j++)
{
c[i, j] = 0;
for (k = 0; k < n1; k++)
{
c[i, j] = c[i, j] + a[i, k] * b[k, j];
}
}
}
result = c;
return result;
}
//二维矩阵乘以一维矩阵
public double[] MatMulti(double[,] a, double[] b)
{
double[] result = null;
int m1 = a.GetLength(0);
int n1 = a.GetLength(1);
int m2 = b.GetLength(0);
if (n1 != m2) return null;
double[] c = new double[m1];
int i, j, k;
for (i = 0; i < m1; i++)
{
c[i] = 0;
for (k = 0; k < n1; k++)
{
c[i] = c[i] + a[i, k] * b[k];
}
}
result = c;
return result;
}
//求矩阵方差
public double MatVaria(double[,] a)
{
int m = a.GetLength(0);
int n = a.GetLength(1);
int i, j;
double sum = 0;
double aeve = MatAvera(a);
for (i = 0; i <= m - 1; i++)
{
for (j = 0; j <= n - 1; j++)
{
sum = sum + Math.Pow((a[i, j] - aeve), 2);
}
}
double result = sum / m / n;
return result;
}
//求方差 int型
public double MatVaria(int[,] a)
{
int m = a.GetLength(0);
int n = a.GetLength(1);
int i, j;
double sum = 0;
double aeve = MatAvera(a);//矩阵均值
for (i = 0; i <= m - 1; i++)
{
for (j = 0; j <= n - 1; j++)
{
sum = sum + Math.Pow((a[i, j] - aeve), 2);
}
}
double result = sum / m / n;
return result;
}
//求矩阵均值
public double MatAvera(double[,] a)
{
double result = 0;
int m = a.GetLength(0);
int n = a.GetLength(1);
int i, j;
double sum = 0;
for (i = 0; i < m; i++)
{
for (j = 0; j < n; j++)
{
sum = sum + a[i, j];
}
}
result = sum / m / n;
return result;
}
//求矩阵均值
public double MatAvera(int[,] a)
{
double result = 0;
int m = a.GetLength(0);
int n = a.GetLength(1);
int i, j;
double sum = 0;
for (i = 0; i < m; i++)
{
for (j = 0; j < n; j++)
{
sum = sum + a[i, j];
}
}
result = sum / m / n;
return result;
}
//求两个矩阵协方差
public double MatCovar(double[,] a, double[,] b)
{
//还需要先判断两个矩阵形式
int m = a.GetLength(0);
int n = a.GetLength(1);
int i, j;
double aev1 = MatAvera(a);
double aev2 = MatAvera(b);
double sum = 0;
for (i = 0; i <= m - 1; i++)
{
for (j = 0; j <= n - 1; j++)
{
sum = sum + (a[i, j] - aev1) * (b[i, j] - aev2);//协方差
}
}
double result = sum / m / n;
return result;
}
public double MatCovar(int[,] a, int[,] b)
{
int m = a.GetLength(0);
int n = a.GetLength(1);
int i, j;
double aev1 = MatAvera(a);
double aev2 = MatAvera(b);
double sum = 0;
for (i = 0; i <= m - 1; i++)
{
for (j = 0; j <= n - 1; j++)
{
sum = sum + (a[i, j] - aev1) * (b[i, j] - aev2);
}
}
double result = sum / m / n;
return result;
}
//求两个矩阵相关系数
public double MatCorr(double[,] a, double[,] b)
{
double var1 = MatVaria(a);
double var2 = MatVaria(b);
double cov12 = MatCovar(a, b);
double result = cov12 / Math.Sqrt(var1 * var2);
return result;
}
public double MatCorr(int[,] a, int[,] b)
{
double var1 = MatVaria(a);
double var2 = MatVaria(b);
double cov12 = MatCovar(a, b);
double result = cov12 / Math.Sqrt(var1 * var2);
return result;
}
//求逆函数
public double[,] MatInver(double[,] n) //矩阵求逆函数 组员法改进
{
//前提判断: 是否为方形矩阵 是否满足可逆条件
int m = n.GetLength(0);
double[,] q = new double[m, m]; //法方程系数逆矩阵;
double u, temp; //临时变量
int i, j, k;
//初始单位阵
for (i = 0; i < m; i++)
for (j = 0; j <= m - 1; j++)
q[i, j] = (i == j) ? 1 : 0;
// 求左下
for (i = 0; i <= m - 2; i++)
{
//提取该行的主对角线元素
u = n[i, i]; //可能为0
if (u == 0)
{
for (i = 0; i < m; i++)
{
k = i;
for (j = i + 1; j < m; j++)
{
if (n[j, i] != 0)
{
k = j;
break;
}
}
if (k != i)
{
for (j = 0; j < m; j++)
{
temp = n[i, j];
n[i, j] = n[k, j];
n[k, j] = temp;
//伴随交换
temp = q[i, j];
q[i, j] = q[k, j];
q[k, j] = temp;
}
}
else
MessageBox.Show("不可逆矩阵", "ERROR", MessageBoxButtons.OK);
}
}
for (j = 0; j < m; j++)//该行除以主对角线元素的值 使主对角线元素为1
{
n[i, j] = n[i, j] / u; //分母不为0
q[i, j] = q[i, j] / u; //伴随矩阵
}
for (k = i + 1; k < m; k++) //下方的每一行减去 该行的倍数
{
u = n[k, i]; //下方的某一行的主对角线元素
for (j = 0; j < m ; j++)
{
n[k, j] = n[k, j] - u * n[i, j]; //下方的每一行减去该行的倍数 使左下角矩阵化为0
q[k, j] = q[k, j] - u * q[i, j]; //左下伴随矩阵
}
}
}
u = n[m - 1, m - 1]; //最后一行最后一个元素
if (u == 0)
MessageBox.Show("不可逆矩阵", "ERROR", MessageBoxButtons.OK);
n[m - 1, m - 1] = 1;
for (j = 0; j < m; j++)
{
q[m - 1, j] = q[m - 1, j] / u;
}
// 求右上
for (i = m - 1; i >= 0; i--)
{
for (k = i - 1; k >= 0; k--)
{
u = n[k, i];
for (j = 0; j < m; j++)
{
n[k, j] = n[k, j] - u * n[i, j];
q[k, j] = q[k, j] - u * q[i, j];
}
}
}
return q;
}
}