How to create randomly points within polygons for each row of a dataframe matching a Polygon ID with R

My question is similar to this post, but because there are no "perfect" answer and above all no reproducible example I open this question.

I have a "SpatialPolygonsDataFrame" and a data frame. Both have a variable of common factors. I would like to generate a random point for each observation of the data frame within the according polygon.

I am starting to learn GIS with R. I read the first chapters of Bivand et al. 2013 but i am still strugling with the basics. I tried with the function spsample from the library(sp) but I am puzzling with the variable n of this function: it's an (approximate) sample size

In my example there are 3 polygons and 30 observations. The common variable is MatchID. I need a SpatialPolygonsDataFrame with 15 points in the first polygons ("abcd"), 10 in the second ("efgh") in 5 in the third ("ijkl"). See

table(df$MatchID)

Example:

 library(sp)
# Definition of the CRS
poly.crs <- CRS("+proj=utm +zone=36 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")


# Definition of 3 Polygons
poly.a <- Polygon((matrix(c(-4.653724,1.2210259,-3.803160, 1.3799638, -3.245480, 0.1665987,-4.666098, 0.1097523, -4.653724, 1.2210259), nrow=5, ncol=2, byrow=T)))
poly.b <- Polygon((matrix(c(-5.820343, 2.675320,-5.427519, 3.062975,-4.701119, 2.819967,-4.555540, 1.855489, -5.050758, 1.388308, -5.783673, 1.572882,-5.820343, 2.675320), nrow=7, ncol=2, byrow=T)))
poly.c <- Polygon((matrix(c(-3.758639, 2.873654, -3.273958, 3.417311, -2.213099, 2.935320, -1.972031, 1.945992, -3.033510, 1.279312, -3.709359, 1.844633, -3.758639, 2.873654), nrow=7, ncol=2,     byrow=T))) 

# Making a SpatialPolygons
polys.a = Polygons(list(poly.a), "pa")
polys.b = Polygons(list(poly.b), "pb")
polys.c = Polygons(list(poly.c), "pc")

Spolys = SpatialPolygons(list(polys.a,polys.b,polys.c), 1:3,  proj4string=poly.crs)

# Making a SpatialPolygonsDataFrame
data.Spolys<- (data.frame(MatchID=c("abcd","efgh","ijkl"), row.names=row.names(Spolys)))
Poly <- SpatialPolygonsDataFrame(Spolys, data.Spolys, match.ID = TRUE)


## Data frame that should be converted in a SpatialPointsDataFrame thanks to the variable "MatchID"
df <- data.frame(
          MatchID=c("abcd","abcd","abcd","efgh","efgh","ijkl"),

          V1 = 1:30,
          V2 = "a"
)

I need a function to generalize the idea below (looking for the matching with the factor) and giving exactly the number of observation I need (and not an approximate random size)

spdf.a <- SpatialPointsDataFrame(spsample(Poly[Poly$MatchID=="abcd",], n = nrow(df[df$MatchID=="abcd",]), "stratified"), df[df$MatchID=="abcd",])

Answer

After a lot of attempts I have this solution, probably not so clean. Comments, improvements or other way to answer are much welcome!

### Preparing the SpatialPointsDataFrame
spdf <- matrix(as.numeric(NA), nlevels(Poly$MatchID), 1)

spdf <- as.list(spdf)

### Sample the coordinate, match it with data in spdf. It creates a list fore each factor of the MatchID
### sample(spsample()) fix the size of the sample
for (i in seq(Poly$MatchID))
    spdf[i] <- SpatialPointsDataFrame(
               sample(spsample(Poly[order(Poly$MatchID)==i,], n = 100, "stratified"),table(df$MatchID)[[i]]),  ### table(df$MatchID)[[i]] is the size of the sample and match the sum of factors in df 
               df[df$MatchID==dimnames(table(df$MatchID))[[1]][i],], ##  dimnames(table(df$MatchID))[[1]][i] ### match the value of the selected "factor" to select the rows of the data
               proj4string=poly.crs, 
               match.ID=TRUE)


## Merging together the list to make a SpatialDataFrame
spdf <- do.call("rbind", spdf)

## Plot 
plot(Poly[,])
plot(spdf, add=TRUE, col=spdf$MatchID)