Implementación de M-spline en R

Question

Estoy implementando M-spline en R como se define aquí : http://www.fon.hum.uva.nl/praat/manual/spline.html y originalmente en Ramsay (1988).

En resumen, definimos una lista de nudos $t$ tal que :

$0=t_1=...=t_k<t_{k+1}<...<t_{k+m}<t_{k+m+1}=...=t_{k+m+k}=1$

Así que el $k$ los primeros nudos son 0 y los $k$ los últimos nudos son 1 con $m$ nudos interiores.

El spline Mi de orden k con nudos t se define recursivamente como sigue :

$M_i(x|1,t) = 1 / (t_{i+1} – t_i), t_i x < t_{i+1},$ 0 en caso contrario

$M_i(x|k,t) = k [(x–t_i)M_i(x|k–1,t) + (t_{i+k}–x)M_{i+1}(x|k–1,t)] / ((k–1)(t_{i+k}–t_i))$

El problema es que si quiero las M-splines de orden 3 (k=3), el $M_1$ spline y el $M_{k+m}$ spline no existen porque cuando se calcula $M_i(x|k–1,t)$ se divide por $(t_{i+k}–t_i)$ pero desde ahora $k=3-1=2$ y $i=1$ (o $m+k$ ), $(t_{i+k}–t_i)$ es $(t_{1+2}–t_1)=0$ ya que los tres primeros nudos son 0 (lo mismo para el último).

Así que la definición parece tener un problema. Aquí está mi implementación en R con un ejemplo del problema. Aquí k=3 y m=3.

ts = c(0,0,0,0.3,0.5,0.6,1,1,1)
Mk = function(i,k,x,ts){
  if(k==1){
    if(ts[i]<=x && x<ts[i+1]){1/(ts[i+1]-ts[i])}
    else{0}
  }
  else{
    #print(paste(i,k))
    #print((ts[i+k]-ts[i]))
    k*((x-ts[i])*Mk(i,k-1,x,ts)+(ts[i+k]-x)*Mk(i+1,k-1,x,ts))/((k-1)*(ts[i+k]-ts[i]))

  }
}
Mk(1,3,.4,ts)

Tenga en cuenta que mi código y la definición parecen funcionar para todas las splines "internas".

Gracias.

Answer 1

Aunque sé que no es una pregunta nueva, sólo quiero mencionar una implementación existente de M-splines en R como referencia.

Paquete splines2 proporciona una función llamada mSpline para las M-splines. Si tiene experiencia en el uso del paquete splines probablemente ha sabido utilizar mSpline ya que su interfaz de usuario es exactamente igual con bs para B-splines en el paquete splines .

Además, las bases de las M-splines se evalúan aprovechando una sencilla transformación entre las bases de las B-splines y las M-splines, mientras que la evaluación de las B-splines se realiza de forma eficiente mediante el paquete splines e implementado en C. Por lo tanto, la función mSpline debería tener un mejor rendimiento en velocidad, en comparación con la implementación directa mediante fórmulas recursivas puramente en R.

Una demostración rápida está disponible en el paquete viñetas .

Answer 2

Bueno, después de algunos retoques con mi código he intentado añadir esta línea a la definición de Mk :

    if(ts[i+k]-ts[i]==0){0}

Para que cuando salga de la lista, la spline sea simplemente cero. Funcionó y pude confirmar que la base era correcta.