CODEWITHSUNDEEP
rmatlabsmoothing
CHP
CHP

Reputation: 950

Smoothing with lowess

I'm trying to convert a matlab script to R, and am having some trouble with smoothing.

The matlab code I want to convert is the following:

for i = 1:size(spike_sum,2)
        smooth_sum(1:Ne,i)=smooth(double(spike_sum(1:Ne,i)),spanNe,'lowess');
end

for i = 1:Ne
        smoother_sum(i,:)=smooth(double(smooth_sum(i,:)),spanT,'lowess');
end

where spike_sum is a matrix that is Ne x 4000. I want to first smooth in Dim 1, with span spanNe, and do that for all 4000 slices. Then, I want to smooth in Dim 2 with span spanT, and do that for all Ne slices.

I've looked at the lowess function in R, but it seems like it takes in 2 dimensions as lowess(x,y,span,iter,delta). So to get the results of the code above in R, would I just take a slice of the matrix for y and replicate a constant value for x?

Upvotes: 1

Views: 2351

Answers (2)

CHP
CHP

Reputation: 950

Well, although the answer provided by joran did not work as intended, and joran spent quite a bit of time trying to debug it, which I appreciate, ultimately I was unable to get that solution to work, and I don't really know why, but it was not producing results comparable to the matlab code. After messing around I found out this solution (which might not be optimal, but works as expected)

loesscontrol=loess.control(surface="interpolate", statistics="approximate", trace.hat="exact", iterations=1)
spanNe=100/Ne
spanT=50/nsteps
spanNi=30/Ni

for (i in 1:nsteps){
        x<-1:Ne
        y<-spike_sum[1:Ne,i]
        smoothingNe<-loess.smooth(x, y, span=spanNe, degree=1, family="gaussian", evaluation=Ne)
        smooth_sumNe[1:Ne,i]<-smoothingNe$y
}

for (i in 1:Ne){
        x<-1:nsteps
        y<-smooth_sumNe[i,1:nsteps]
        smoothingT<-loess.smooth(x, y, span=spanT, degree=1, control=loesscontrol, family="gaussian", evaluation=nsteps)
        smoother_sumNe[i,1:nsteps]<-smoothingT$y        
}

I want to mention that one of the keys here was to set evaluation=Ne in the first smoothing, because otherwise the result was null. I don't know why that was, but perhaps because the data was sparse and highly discontinuous.

Upvotes: 0

joran
joran

Reputation: 173627

My Matlab is quite rusty, but if I've understood you correctly, you probably want to pass the sequence 1:Ne or 1:4000 for the x argument to lowess, as that is meant to be the x coordinates of the points your are smoothing. This assumes that you are assuming your points are indeed equally spaced.

Something like this would work:

#Example matrix
M <- matrix(runif(1600),40,40)

#Smooth rows; transpose when smoothing over rows
M1 <- t(apply(M,1,FUN = function(x){lowess(1:length(x),x)$y}))

#Smooth columns; but don't transpose; fills by column already
M2 <- apply(M1,2,FUN = function(x){lowess(1:length(x),x)$y})

I didn't include the different spans, but you can add those details in yourself.

Edit

Ah, but if you're looking for speed, you should probably be using loess.smooth directly. loess uses a formula interface, so you'll want to call loess.smooth directly. It's defaults are different than lowess, though, so be careful. Swapping that function in cut the running time for me by almost 1/4.

Upvotes: 1

Related Questions