Project 1 of Statistics.md

一、准备阶段

本次实验主要考虑影响大倪体重各体型指标的因素。我们的数据收集了112组大鲵的体重(weight)、全长(tl)、体长(bl)、头长(hl)、体高(h)、体宽(w)、尾柄长(cl)和肠长(sl)。

利用统计软件R来实现,其中用到的包有readxl, dplyr, ggplot2, car

阅读更多
Assignment 1 of Advanced Computational Method

Homework 1

Monte Carlo method is a method to compute the expectation of a random variable.

The basic steps follow like this:

  1. draw random samples ${x_1, x_2, \cdots, x_n}$ according to $p(x)$.

  2. derive the monte carlo estimator of $E(f)$ from

$$
\hat{f}{\text{mc}}=\frac{1}{n}\sum{i=1}^n f(x_i)
$$

阅读更多
Numerical Examples in Java

Some algorithms from the textbook Numerical Analysis by Richard L. Burden and J.Douglas Fairs (9th ed)..

DividedDifference.java

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
package Langrange;

import java.util.Arrays;

public class DividedDifference {

private static double dividedDiffHelp(double[] x, double[] y, int n, int i) {
//compute the nth divided difference
if (n == 0) {
return (y[i]);
} else {
return (dividedDiffHelp(x, y, n - 1, i + 1) - dividedDiffHelp(x, y, n - 1, i)) / (x[i + n] - x[i]);
}
}

public static double[] dividedDiff(double[] x, double[] y) {
//compute divided diff needed for each coefficient
double[] coef = new double[x.length];
for (int i = 0; i < x.length; i++) {
coef[i] = dividedDiffHelp(x, y, i, 0);
}
return coef;
}

public static String dividedOut(double[] x, double[] y) {
String out = "";
double[] coef = dividedDiff(x, y);
for (int i = 0; i < coef.length; i++) {
String temp = "";
if (i == 0) {
temp = "" + coef[0];
} else {

for (int j = 0; j < i; j++) {
temp = temp + "(x-" + x[j] + ")";
}
temp = coef[i] + temp;
}
if (i != coef.length - 1) {
temp = temp + "+";
}
out = out + temp;

}
return out;
}

public static String printTable(double[] x, double[] y) {
//print line by line
// number of columns will be the number of points -1+2
String out = "";
for (int row = 0; row < x.length; row++) {
out = out + row + " " + x[row] + "\t" + y[row] + "\t";
for (int col = 1; col <= row; col++) {
out = out + dividedDiffHelp(x, y, col, row - col) + "\t";
}
//new row
out = out + "\n";
}
return out;
}

public static double interpolate(double x, double[] coef, double[] data) {
double y = 0;
for (int i = 0; i < coef.length; i++) {
double product = 1;
for (int j = 0; j < i; j++) {
product = product * (x - data[j]);
}
y = y + product * coef[i];
}
return y;
}

public static void main(String[] args) {
//PROBLEM 8
//test dividedDiff
double[] x = { 0, 0.1, 0.3, 0.6, 1.0, 1.1 };
double[] y = { -6, -5.89483, -5.65014, -5.17788, -4.28172, -3.99583 };
//System.out.println(dividedDiffHelp(x,y,4,0));
// double[] coef = dividedDiff(x,y);
//System.out.println(Arrays.toString(coef));
System.out.println(dividedOut(x, y));

//PROBLEM 9
double[] x9 = { 0, 0.4, 0.7 };
double[] y9 = { 1, 3, 6 };
System.out.println(dividedOut(x9, y9));
System.out.println(printTable(x9, y9));
System.out.println(interpolate(0, dividedDiff(x9, y9), x9));
}

}

NaturalCubicSpline.java

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
package CubicSpline;

import java.util.Arrays;

public class NaturalCubicSpline {
//this class file will compute the natural cubic spline
//we are interpolating over data array x with values y

public static double[][] fitSpline(double[] x, double[] y) {
//compute hi=xi+1-xi
double[] h = new double[x.length - 1];
for (int i = 0; i < h.length; i++) {
h[i] = x[i + 1] - x[i];
}
double[] a = new double[h.length];
//a[0] left empty
for (int i = 1; i < a.length; i++) {
a[i] = 3 / h[i] * (y[i + 1] - y[i]) - 3 / h[i - 1] * (y[i] - y[i - 1]);
}
//create extra arrays
double[] l = new double[x.length];
double[] u = new double[x.length];
double[] z = new double[x.length];
l[0] = 0;
z[0] = 0;
for (int i = 1; i < x.length - 1; i++) {
l[i] = 2 * (x[i + 1] - x[i - 1]) - h[i - 1] * u[i - 1];
u[i] = h[i] / l[i];
z[i] = (a[i] - h[i - 1] * z[i - 1]) / l[i];
}
l[l.length - 1] = 1;
z[z.length - 1] = 0;
double[] c = new double[x.length];
c[c.length - 1] = 0;

double[] b = new double[x.length - 1];
double[] d = new double[x.length - 1];

for (int j = x.length - 2; j >= 0; j--) {
c[j] = z[j] - u[j] * c[j + 1];
b[j] = (y[j + 1] - y[j]) / h[j] - h[j] * (c[j + 1] + 2 * c[j]) / 3;
d[j] = (c[j + 1] - c[j]) / 3 / h[j];
}

double[][] coef = new double[x.length - 1][4];
for (int j = 0; j < x.length - 1; j++) {
coef[j][0] = y[j];
coef[j][1] = b[j];
coef[j][2] = c[j];
coef[j][3] = d[j];
}

return coef;

}

public static double interpolate(double x, double[] coef, double xdat) {
double val = 0;
for (int i = 0; i < coef.length; i++) {
val = val + coef[i] * Math.pow((x - xdat), i);
}
return val;
}

public static void main(String[] args) {
double[] x = { 1, 2, 3 };
double[] y = { 2, 3, 5 };
double[][] coef = fitSpline(x, y);
System.out.println(Arrays.toString(coef[1]));
System.out.println(interpolate((double) 2, coef[0], x[0]));

//Assignment numerical 7
double[] years = { 1950, 1960, 1970, 1980, 1990, 2000 };
double[] pop = { 151326, 179323, 203302, 226542, 249633, 281422 };
double[][] coef2 = fitSpline(years, pop);
System.out.println(Arrays.toString(coef2[0]));
System.out.println(interpolate(1940.0, coef2[0], years[0]));
System.out.println(Arrays.toString(coef2[2]));
System.out.println(interpolate(1975, coef2[2], years[2]));
System.out.println(Arrays.toString(coef2[4]));
System.out.println(interpolate(2020, coef2[4], years[4]));
}

}