请高手指点以下t分布最大似然估计的算法：

function va=ml_t(data)
    [m,n]=size(data);
    n=max(m,n);
    syms v x;
    tmp = gamma((v+1)/2)/gamma(v/2);
    f = tmp/(sqrt(v*pi)*(1+(x^2)/v)^((v+1)/2));
    log_f=log(f);
    dlog_f=diff(log_f,'v');
    L=sym('0');
    for k=1:n,
        l=subs(dlog_f,x,data(k));
        L=L+l;
    end
    va=solve(L);

  这是我自己编的程序。我用trnd命令生成数据验证过，例如用trnd(3,30,1)生成30个服从自由度为3的t分布的数据，再利用这个函数估计参数是否接近3。结果是这套程序能够顺利通过。但是遇到以下问题：

1、程序运行过程非常非常慢！我逐条语句检查过，“瓶颈”在for循环和最后的solve这两步，而这两处之所以出现“瓶颈”是因为“l=subs(dlog_f,x,data(k));”的结果本身就相当复杂，经过几十次循环相加之后就更加复杂，这又导致最终的L相当复杂，解起来很费劲。

2、结果似乎大大地偏离了原来设定的自由度。例如上面提到的用trnd(3,30,1)生成30个服从自由度为3的t分布的数据，结果似乎远离3！（不过由于每次运行这个程序都需要花很多时间，试得不多，而且每次的数据量都不大）

我希望能够简化这套程序，不知可行不可行？或者有没有其他更好的思路？