randu dataset, part 2

Bogumił Kamiński

10 years ago

[This article was first published on R snippets, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

In my last post I have plotted randu dataset to show that all its points lie on 15 parallel planes. But I was not fully satified with the solution and decided to show this numerically.

It can be done in four steps:

identifying four points lying on the same plane and finding its equation (we know that we have 15 planes so it is enough to check 15*3+1 points to find them);
applying the equation to all points in randu dataset to divide them into 15 classes;
veryfying that the classes separate points into 15 parallel planes;
plotting the solution changing colors for points on different planes.

Here is the code I used to do this:

library(rgl)< o:p>

library(caTools)< o:p>

# STEP 1< o:p>

all.combs <- combs(1:46, 4)< o:p>

i <- 1< o:p>

repeat {< o:p>

model <- lm(z ~ x + y, data = randu[all.combs[i, ], ])< o:p>

if (summary(model)$r.squared > 0.99999) {< o:p>

break< o:p>

}< o:p>

i <- i + 1< o:p>

}< o:p>

# STEP 2< o:p>

line.class <- predict(model, randu) – randu$z< o:p>

line.class <- factor(round(line.class) + 10)< o:p>

# STEP 3< o:p>

summary(lm(z ~ x + y + line.class, data = randu))< o:p>

# STEP 4< o:p>

with(randu, plot3d(x, y, z, axes = FALSE, col = line.class,< o:p>

xlab = “”, ylab = “”, zlab = “”))< o:p>

rgl.viewpoint(theta = –3.8, phi = 3.8, fov = 0, zoom = 0.7)< o:p>

At step three we can see that the model obtains perfect fit. The final figure is plotted below:

To leave a comment for the author, please follow the link and comment on their blog: R snippets.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.