在神经科学研究领域,数据的处理和分析是至关重要的。随着技术的进步,现在可以使用.NET语言来构建中间件应用程序,以处理和预测来自脑机接口(BCI)游戏应用的功能性磁共振成像(fMRI)和脑电图(EEG)数据。本文是系列文章的第一部分,探讨了使用开源库构建基于.NET的时间序列应用程序,这些应用程序可以用于神经科学研究以及商品/股票价格预测。目标是提供开源的非线性时间序列算法作为统计引擎的一部分,例如TIME JAQUAR,用于时域和频域,能够在多CPU/多GPU平台上运行,以处理来自BCI设备的实时建模和预测。
本文不涉及时间序列分析的教程,因为网络上已经有很多优秀的教程,以及作者之前所写的内容。本文主要介绍在.NET中用于原型设计的一些优秀库,这些库可以作为命名空间和原型模板,帮助开发线性和非线性时间序列应用程序。
以下是在应用程序中使用的两个库:
有关这些库的更多信息,请查阅相关文档。当然,线性代数和优化方法的可扩展性是基础,也是后续文章的主题。作者在图1中提供了基础信息,以说明两个库的一些特性,并将代码发布在CodePlex上,以便读者可以跟进项目。
以下是进行估计的步骤:
图1. 步骤1-4的代码列表:
C#
using ABMath.ModelFramework.Models;
using ABMath.ModelFramework.Data;
using MathNet.Numerics.Distributions;
using MathNet.Numerics.RandomSources;
using MathNet.Numerics;
namespace TCW_ModelMaker3
{
class Program
{
static void Main(string[] args)
{
runARMATest();
}
public static void ConsoleWrite()
{
Console.WriteLine("\nWelcome to ModelMaker 3");
Console.ReadKey();
}
public static void runARMATest()
{
// Specify the AR and MA parameters
ARMAModel model1 = new ARMAModel(4, 3);
ARMAModel model2 = new ARMAModel(0, 1);
// Select an random distribution for the data
var sd = new StandardDistribution();
// Make new time series variable
TimeSeries ts1 = new TimeSeries();
var dt = new DateTime(2001, 1, 1);
var normalSim = new StandardDistribution();
var current = new DateTime(2001, 1, 1);
// Create the data
for (int t = 0; t < 100; ++t)
{
double s1 = normalSim.NextDouble();
ts1.Add(current, s1, false);
current = new DateTime(current.AddDays(1).Ticks);
}
model1.SetInput(0, ts1, null);
// Maximum Likelihood Estimation
model1.FitByMLE(10, 10, 10, null);
// Compute the residuals
model1.ComputeResidualsAndOutputs();
// Get the predicted values
var preds = model1.GetOutput(3) as TimeSeries;
var predName = model1.GetOutputName(3);
// Write the model description with parameter values
Console.WriteLine(model1.Description);
Console.ReadKey();
}
}
}
上述代码中最有趣的方法是FitByMLE方法。这个函数使用Halton序列从参数空间中采样,并通过对数似然优化来获得最佳模型。参数值有三种状态:ParameterState.Locked、ParameterState.Free或ParameterState.Consequential。锁定参数在优化方法中保持当前值,自由参数被优化,而结果参数则作为其他参数和数据的函数计算。有关更多信息,请查阅Cronos库。
C#
public virtual void FitByMLE(int numIterationsLDS, int numIterationsOpt, double consistencyPenalty, Optimizer.OptimizationCallback optCallback)
{
thisAsMLEEstimable = this as IMLEEstimable;
if (thisAsMLEEstimable == null)
throw new ApplicationException("MLE not supported for this model.");
int optDimension = NumParametersOfType(ParameterState.Free);
int numConsequential = NumParametersOfType(ParameterState.Consequential);
int numIterations = numIterationsLDS + numIterationsOpt;
var trialParameterList = new Vector[numIterationsLDS];
var trialCubeList = new Vector[numIterationsLDS];
var hsequence = new HaltonSequence(optDimension);
for (int i = 0; i < numIterationsLDS; ++i)
{
Vector smallCube = hsequence.GetNext();
Vector cube = CubeInsert(smallCube);
trialParameterList[i] = thisAsMLEEstimable.CubeToParameter(cube);
trialCubeList[i] = cube;
}
var logLikes = new double[numIterationsLDS];
for (int i = 0; i < numIterationsLDS; ++i)
{
Vector tparms = trialParameterList[i];
if (numConsequential > 0)
{
tparms = ComputeConsequentialParameters(tparms);
trialParameterList[i] = tparms;
}
double ll = LogLikelihood(tparms, consistencyPenalty);
logLikes[i] = ll;
if (optCallback != null)
lock (logLikes)
optCallback(tparms, ll, i * 100 / numIterations, false);
}
// Step 1: Just take the best value.
Array.Sort(logLikes, trialParameterList);
Parameters = trialParameterList[numIterationsLDS - 1];
// Step 2: Take some of the top values and use them to create a simplex, then optimize
// further in natural parameter space with the Nelder Mead algorithm.
// Here we optimize in cube space, reflecting the cube when necessary to make parameters valid.
var simplex = new List();
for (int i = 0; i <= optDimension; ++i)
simplex.Add(FreeParameters(thisAsMLEEstimable.ParameterToCube(trialParameterList[numIterationsLDS - 1 - i])));
var nmOptimizer = new NelderMead { Callback = optCallback, StartIteration = numIterationsLDS };
currentPenalty = consistencyPenalty;
nmOptimizer.Minimize(NegativeLogLikelihood, simplex, numIterationsOpt);
Parameters = ComputeConsequentialParameters(thisAsMLEEstimable.CubeToParameter(CubeFix(CubeInsert(nmOptimizer.ArgMin))));
ComputeResidualsAndOutputs();
}
当然,这是基于Nelder Mead的主要算法,也是未来文章的主题。但这应该足以让开始。