Provided I understood your question correctly, here is a suggestion. I made up some data for population, so if those are different for you the aggregate function call might need to be adapted.
library(rgdal)
library(sp)
# shapefile
tmpdir <- tempdir()
download.file("http://www.abs.gov.au/ausstats/subscriber.nsf/log?openagent&1259030001_ste11aaust_shape.zip&1259.0.30.001&Data%20Cubes&D39E28B23F39F498CA2578CC00120E25&0&July%202011&14.07.2011&Latest", "ste11aaust.zip")
unzip("ste11aaust.zip", exdir = tmpdir )
ste <- paste0(tmpdir,"/STE11aAust.shp")
AUS <- readOGR(ste, "STE11aAust")
unlink(tmpdir)
# make up the population data
set.seed(100)
lon <- runif(10, 1, 147)
lat <- runif(10, -31, -20)
no_ppl <- sample(1:100, 10)
# create SpatialPointsDataframe with people
people.spdf <- SpatialPointsDataFrame(cbind(lon, lat), data.frame(no_ppl), proj4string=CRS(proj4string(AUS)))
# find which polygons the points fall into and sum them up
ppl.sum <- aggregate(x = people.spdf["no_ppl"], by = AUS, FUN = sum)
# plot crude choropleth map (note that there are NAs for some states)
spplot(ppl.sum, "no_ppl")
