using UnityEngine;
using MathNet.Numerics;
using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.Optimization;
using System;
public class MathNETController : MonoBehaviour
{
public GameObject PointsPrefab;
GameObject line;
Vector<double> x;
Vector<double> y;
float scale;
int n;
int count;
// Start is called once before the first execution of Update after the MonoBehaviour is created
void Start()
{
var obj = ObjectiveFunction.NonlinearModel(OptimizeFunction, Gradient, Xdata, Ydata);
var solver = new LevenbergMarquardtMinimizer();
var result = solver.FindMinimum(obj, InitialValue1);
x = Vector<double>.Build.DenseOfArray(Generate.LinearSpaced(1000, 1.0, 15.0));
y = OptimizeFunction(result.MinimizingPoint, x);
line = GameObject.Find("Line");
scale = 50.0f;
n = x.Count;
count = 0;
for (int i = 0; i < Xdata.Count; i++)
{
GameObject point = Instantiate(PointsPrefab);
point.transform.position = new Vector3((float)Xdata[i], Convert.ToSingle(Ydata[i]/scale), 0f);
}
}
// Update is called once per frame
void Update()
{
if (count < n) {
line.transform.position = new Vector3((float)x[count], Convert.ToSingle(y[count] / scale), 0f);
count++;
}
}
Vector<double> Xdata = Vector<double>.Build.DenseOfArray(new double[] {
1.00, 2.00, 3.00, 4.00, 5.00, 6.00, 7.00, 8.00, 9.00, 10.00,
11.00, 12.00, 13.00, 14.00, 15.00
});
Vector<double> Ydata = Vector<double>.Build.DenseOfArray(new double[] {
16.08, 33.83, 65.80, 97.20, 191.55, 326.20, 386.87, 520.53, 590.03, 651.92,
724.93, 699.56, 689.96, 637.56, 717.41
});
Vector<double> InitialValue1 = Vector<double>.Build.DenseOfArray(new double[] { 100, 10, 1, 1 });
Vector<double> OptimizeFunction(Vector<double> p, Vector<double> x)
{
var y = Vector<double>.Build.Dense(x.Count);
for (int i = 0; i < x.Count; i++)
{
y[i] = p[0] / Math.Pow(1.0 + Math.Exp(p[1] - p[2] * x[i]), 1.0 / p[3]);
}
return y;
}
private Matrix<double> Gradient(Vector<double> p, Vector<double> x)
{
var prime = Matrix<double>.Build.Dense(x.Count, p.Count);
for (int i = 0; i < x.Count; i++)
{
prime[i, 0] = 1.0 / Math.Pow(1.0 + Math.Exp(p[1] - p[2] * x[i]), 1.0 / p[3]);
prime[i, 1] = -p[0] * Math.Exp(p[1] - p[2] * x[i]) /
(p[3] * Math.Pow(1.0 + Math.Exp(p[1] - p[2] * x[i]), 1.0 / p[3] + 1.0));
prime[i, 2] = p[0] * x[i] * Math.Exp(p[1] - p[2] * x[i]) /
(p[3] * Math.Pow(1.0 + Math.Exp(p[1] - p[2] * x[i]), 1.0 / p[3] + 1.0));
prime[i, 3] = Math.Log(1.0 + Math.Exp(p[1] - p[2] * x[i])) * p[0] /
(p[3] * p[3] * Math.Pow(1.0 + Math.Exp(p[1] - p[2] * x[i]), 1.0 / p[3]));
}
return prime;
}
}
最近のコメント